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Using Fourier Based Modulating Signals and High Order Statistics by JianQiang 
Pan, Ph.D., Brown University, May 1992. 


This research investigates several important problems in the fields ol signal pro- 
cessing and model identification, such as system structure identification, frequency 
response determination, high order model reduction, high resolution frequency anal- 
ysis, deconvolution filtering, and etc. Each of these topics involves a wide range of 
applications and has received considerable attention. 

Using the Fourier based sinusoidal modulating signals, it is demonstrated that a 
discrete autoregressive model can be constructed for the least squares identification 
of continuous systems. Some identification algorithms are presented lor both SISO 
and MIMO systems frequency response determination utilizing only transient data. 
Also, several new schemes for model reduction have been developed. 

Based upon the complex sinusoidal modulating signals, a paramctiic least squares 
algorithm for high resolution frequency estimation is proposed. Numerical exam- 
ples demonstrate that the proposed algorithm shows better performance than the 
well-known High-Order Yule- Walker Estimation. Also, we have studied the problem 
of deconvolution and parameter identificaton of a general noncausal nuimiinimum 
phase ARMA system driven by non-Gaussian stationary random processes. Inverse 
cumulant and inverse polyspectra are introduced as generalizations of the inverse au- 
tocorrelation and inverse power spectrum. Algorithms are presented for estimating 
the inverse cumulants, both in the frequency domain via the FFT algorithms and in 
the time domain via the least squares algorithm. Simulation results are provided to 
demonstrate the performance of the method. 
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Chapter 1 


Introduction 


This research investigates several important problems in the fields of signal filtering, 
model identification, such as system and structure identification, frequency response 
determination, high-order model reduction, high resolution frequency analysis, de- 
convolution filtering, and etc. All these topics involve a wide range of applications 
and have received considerable attention. 

System modeling and structure identification play a very important role in process 
control and analysis. Much investigation has been carried out in modeling linear 
systems, both in a deterministic and stochastic vein. In a deterministic setting, we 
can use the classical state steady frequency domain approaches , or the time domain 
approaches, like the Bellman-Kalaba quasilinearization technique [Bell65, Kala82], 
model adaptive observer , state variable filters , and etc. In a stochastic setting, vve 
could use the well-known Kalman filter technique , instrumental variable method, 
least squares prediction error approach , and etc.[Astr81, YounSl]. 


1 


In chapter 2, the modulating function approach is formulated for system iden- 
tification, including order determination. Using the Fourier based commensurable 
sinusoids as modulating functions, it is demonstrated that a discrete autoregressive 
model for the least squares identification of continuous differential systems can be con- 
structed in a way utilizing the computationally efficient FFT technique while avoiding 
the complicated problem of estimating unknown initial values of the system. We dis- 
cuss the issue of input persistent excitation based upon the autoregressive model. It 
is shown how an AR model can be constructed to estimate the system parameters 
and to determine the order of the system. It is proposed that the choice of sinusoidal 
modulating functions can be based on the heuristic notion of system bandwidth . 

Determination of the frequency response of a stable linear system from input- 
output data is a classical problem of signal processing and system control, and is very 
important for system synthesis, controller design and the Nyquist stability analysis. 
Methods for solving this problem include the DFT technique, the cross correlation, 
and etc.[Astr75, Ljun87, Unbe87]. Requirements for employing these nonparametric 
schemes are the statistical stationarity of the input and output data, or the steady 
state operating condition of the system. In order to satisfy these assumptions, elim- 
inate the initial value effects and achieve good accuracy, a long data record needs 
to be collected. Good reviews of these approaches, as well as noise and finite data 
effects, can be found in [Sode89] . 

A frequent critical situation is that the model is not initially at rest and the 
available length of input and output signals is limited. Also the signal involves sensor 
noise. In this case, the DFT and cross correlation algorithms will entail large error. 
Chapter 3 deals with this more general situation by using the complex sinusoidal 
modulating signals, see [PearAE, Pear91]. An algorithm for MIMO systems frequency 
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response determination has been developed. Numerical examples are also presented 
to show the performance. 

Based upon the algorithms developed in chapter 3, chapter 4 proposes an ap- 
proach for model reduction. The new method combines the model identification and 
reduction together. There is no necessity of knowing the original system or assuming 
statistical stationarity of input-output data which is crucial to existing algorithms. 
Some comparision with the balanced realization method and the P-L method have 
been made. 

Chapter 5 focuses on the problem of high resolution frequency estimation by us- 
ing short time data and a simple linear least squares algorithm. An autoregressive 
differential equation model is constructed to fit the received signal. Using the Fourier 
based complex sinusiodal modulating signals, we can easily convert the differential 
equation model into linear algebraic equations. Numerical examples demonstrate that 
the proposed algorithm shows better performance than the well-known Fiigh-Order 
Yule- Walker Estimation. 

In chapter 6, we focus on the problem of deconvolution and parameter identificaton 
of a general noncausal nonminimum phase ARMA system driven by non-Gaussian 
stationary random processes. Inverse cumulant and inverse polyspectra are intro- 
duced as generalizations of the inverse autocorrelation and inverse power spectrum 
proposed in [Chat79, Clev72]. The original noncausal ARMA system is approximated 
by a noncausal AR model; a direct relationship exists between the inverse cumulants 
and the parameters of the deconvolution AR filter. The inverse system is then con- 
structed directly by utilizing a gradient type nonlinear optimization algorithm for 
matching the AR coefficients and the estimated inverse cumulants. Algorithms are 
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presented for estimating the inverse cumulants, both in the frequency domain via the 
FFT algorithms and in the time domain via the least squares algorithm. In the time 
domain estimation, we have used the forward regression orthogonal algorithm pro- 
posed in [Bill89] because it can be modified to take advantage of the specific problem 
structure. Simulation results are provided to demonstrate the performance of the 
method. 
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Chapter 2 


Input Persistent Excitation and Model 
Structure Estimation 


2.1 Linear Systems and Moment Functionals 


In this section, we focus on linear systems and introduce the modulating functionals 
which can be used for system identification. A general input-output linear differential 
system can be specified by : 

2 >„-.yy (0 = K-i-ip'u(t) (2.1.1) 

«=0 t =0 

The object is to estimate the structure and parameters of the system given input- 
output data (u(t),y(t)) on the time interval [0, T). 

Shinbrot’s method of moment functionals is one of the classical techniques for sys- 
tem identification[Shin57]. It is also called the modulating function approach. We use 
modulating functions to convert a differential equation into algebraic equations which 
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makes it easier to solve the identification problems. Shinbrot’s modulating function 
technique avoids dealing with the unknown initial conditions over each time sequence 
interval [f,,t, + i], and avoids differentiating the original data , thereby avoiding the 
noise sensitivity in estimating time derivatives of data. 

Converting the differentia] equation into an algebraic equation is a big step, the 
latter is much easier to handle. As introduced by Shinbrot, <j>(t) is a modulating 
function of order N over a fixed time interval [0,T] if 4>{t) is sufficiently smooth and 
satisfies the end point conditions : 

«p (,) (0) = 4>('\T) = 0, * = 0, 1, . . . , N - 1 

where, 

<£ (,) (0 = P'4>{t) = <£(<) 


The significance of using this property for the model identification relates to the 
following fact. If equ.(2.1.1) is multiplied by <j>{t) on both sides and integrated by parts 
n times over [0, T], while using the end point conditions, the result is the following 
functional equation: 


in-lj’an-i [ T <p {,) (t)y(t)dt = "f (-1 YK-i />(<)«(*)* (2.1.2) 

f^o Jo f^o J o 


Here we assume is a N th order modulating function, N > n. Furthermore, if 
i = 1,2,..., 2 M + 1 is a set of linearly independent modulating functions, 
then a set of algebraic equations result . From these vector algebraic equations we can 
use standard least squares techniques to estimate the parameters. 


Several researchers have already investigated similar topics using a variety of mod- 
ulating functions such as Hermite polynomials and splines. However, using these 
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modulating functions has not significantly reduced the computational cost because 
thev lack some kind of fast algorithm. Motivated by this fact we use modulating 
functions comprised of linear combinations of commensurable sinusoids because the 
FFT algorithm can be applied to calculate the following integrals efficiently: 

rT 27T 

J Z(f )[sin(mu>ot) or cos(mojot)]dt , u>o = — 

especially for large M. Note that u> 0 plays the role of a “resolving frequency ” in the 
system idenfication scheme. The saving in computation time is significant, and the 
use of the FFT makes the procedure possibly to be an on-line algorithm. 

With this advantage of using the FFT fast algorithm, we consider the set of 
commensurable sinusoids {cos(mwof), sin(mu;of)} defined on the time interval [0, T], 
u Q = 2x/T, m = 0,1,. . . , M. Let /(f) be the sinusoid vector defined as: 

/(f) = col[l, cos w 0 f, sin u> 0 f,---, cos Muj 0 t, sin Mu> 0 t] (2.1.3) 

It could be expected that if we form the appropriate linear combinations of these 
sinusoids subject to the end point conditions then we can get a set of linear indepen- 
dent trigonometric modulating functions. 

For constructing the N th order modulating functions, using the off-line procedure 
specified in [LeeF84], we can get 2M+l-iV linearly independent modulating functions 
i = 1,2,..., 2 M + 1 — N, which can be represented as the vector-matrix form: 

*(<) = Cf(t), 0 <t<T (2.1.4) 

where C is a (2 M + 1 - N) by (2 M + 1) matrix constructed to make <t>(f) satisfy the 
end point conditions: 

*W(0) = * (i \T) = 0, z = 0,l,...,JV-l (2.1.5) 
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The rows of matrix C are determined by the solution of a Vandermonde type matrix 
equation. Moreover, the matrix C has full rank . 


Because is a vector of the N th order modulating function, from equations 
(2.1.3) and (2.1.4), the time derivatives of $(£) can be easily computed in the form 
of : 

$«(i) = (-ipm * = 0,1,2,... (2.1.6) 

where D is a matrix defined by the block diagonal structure: 

/0 \ 

0 1 

-1 0 


D = 


U/Q 


V 

We assume D° = /, the identity matrix. 


0 M 
—M 0 J 


(2.1.7) 


For identifying a linear time invariant system (2.1.2), multipling both sides by $(<), 
and integrating by parts n times results in the following algebraic vector equation: 

rT ... ,l T_ 1 rT 


D-l = E(-l I 

«'= 0 «=0 


a o = l 


(2.1.8) 


Noting equ.(2.1.6), (2.1.8) can be expressed as : 


n— 1 


Y,a n -iCD'Y = £ K-x-iC&U 


1=0 


1=0 
Clo — 1 


(2.1.9) 
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( U,Y ) are corresponding finite Fourier series coefficient vectors of the input and 
output signals: 

U = I T u{t)f(t)dt 

Jo 

Y = [ T y(t)f(t)dt (2.1.10) 

JO 

Define a parameter vector 6 and coefficient matrix r as : 

6 = co/[-ai, . . . , — a n , 6 0 , • • • ,K-i] 

r = [D n - 1 Y,...,Y,D n - 1 U,...,U] (2.1.11) 

Then, equ.(2.1.9) can be written equivalently in a regression form as : 

CvB = CD n Y (2.1.12) 

We can use the standard least squares technique to solve equation (2.1.12). The 
dimension of the parameter vector 0 for (2.1.12) is 2n. The normal equation is: 

(CT) t (CT)6> = -(Ct) t C D n Y (2.1.13) 


If we choose (Ai,N) such that : 

2M + l-N>2n (2.1.14) 

then there exists a unique solution for a one-shot least squares estimation if and only 
if we can find input-output data for which CT has full rank. Although the matrix C 
has full rank 2 M + 1 - TV, it does not guarantee that CT will have full rank. Thus 
we have to consider C and r together. 

Now there are still some questions which have not been answered. For any general 
system, can we always find some conditions which involve input signals only and 
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guarantee that CT has full rank? Can we find some good input signal which enable 
us to use some efficient algorithms to estimate 6 in a robust way? How can we use 
equ.(2.1.12) to estimate the order of a system? All these questions are related to the 
existence and uniqueness of a one-shot estimation. These issues will obviously involve 
the notions of system identifibility, persistent excitation of the system active modes, 
appropriate structure of the model, and etc. In the following sections we mainly 
investigate these related issues. 


2.2 System Persistent Excitation 


If we choose the Fourier based modulating functions (2.1.4), then equ.(2.1.2) becomes 
(2.1.12) and there exists a unique least squares estimation iff the matrix CT has full 
rank. From a practical point of view, we are interested in some conditions which 
involve the input u(t) only. 

The idea is quite simple. If 4>(f) is an N th order modulating function, then & k \t) 
is also an (N — k) th order modulating function for any k, N — k > n. Based on 
the off-line algorithm [Pear85], for fixed ( M,N ), we form the N ih order modulating 
functions by computing the (2 M + 1 - N) x (2 M + 1) matrix C in equ.(2.1.4). For 
the linear time invariant system(2.1.1), multipling both sides by <*>W(£), and taking 
intergration-by-parts n times, the results are the algebraic equations: 

j2a n -iCD i+k Y = £ CD' +k U 

i=0 i = 0 

<z 0 = 1 , fc = 0, 1 , . . . , JV — n (2.2.1) 

where D,Y,U are defined in the equ.(2.1.7) and (2.1.10). 
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Define: 


W(m) = CD m Y, 
V(m) = C D m U 1 


( 2 . 2 . 2 ) 


then (2.2.1) can be written as: 

(f>,<r)W(rn) = (E ~ 1) 


n— 1 


i=0 i=0 

m = n, n , + 1, . . . , N 


(2.2.3) 


where q~ l is the unit delay operator, i.e. q l W{m) = W(m - 1). 


In equ. (2.2.3) W(m) depends on the output signal, V(m) depends on the input 
signal. The dimension of both W(m) and V(m) is 2 M + 1 — N. So instead of 
using the n th order modulating functions and getting one set of algebraic equations 
having high dimension of order 2A1 + 1 — n, here we use the order modulating 
functions and obtain a sequence of algebraic equations but having the low dimension of 
2 M + 1 — N . Furthermore, if we choose N = 2 M, we get a sequence of scalar algebraic 
equations. The essential information about the system parameters contained in these 
two algebraic equation sets is the same . The advantage of using (2.2.3) is that it 
has the sequence structure that makes it easy to get some conditions on good input 
signals and develop some recursive algorithms for identifing the system parameters 
and estimating the order of the system in a robust manner. 

Realizing that if the data involve some measurement noises, equ. (2.2. 3) should be 
replaced by the linear regression model: 

W(m) = ip T (m)0 + e(m) (2.2.4) 

where: 

0 — (uj, i Ct ^1 j • * * ’ b n — \ ) 
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0(m) = ( — W(m — 1), —W{m — 2), . . . , - W(m- n), V^ro — 1), . . . , V{m — n)) T (2.2.5) 


Denote: 

Mq~ l ) = Y a ^~' 

i=0 

Sfa -1 ) = (2.2.6) 

1=0 

Now, the standard model persistent excitation analysis tool can be used for the 
linear regression model (2.2.4). See [Sode89] for more details regarding the analysis 
we employ in the following. 

The condition that there exists a unique solution to the least squares estimation 
problem of (2.2.4) will be asymptotically equivalent to the nonsingularity of a conva- 
riance matrix: 

R = E^4> t > 0 (2.2.7) 

i.e., R must be positive definite. So R > 0 will be relevant to the persistent excitation 
property of input. Actually, R can be estimated using the following equation: 

R = Y, 4>{m)4> T {m) 

n m=n 

Now we consider two different cases, e(m) = 0, i.e., the system is free of noise, and 
e(m) ^ 0, i.e., the system is contaminated by a significant measurement noise. 

Case 1. e(m) = 0 
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For t = 0, we can make the following analysis similar to [Sode89]. Rewrite: 

f-h ... -b n - , \ 


i/>(m) = 


b\ ... ^n— 1 


/ ]4^=TjV(m - !) \ 




1 ... Cln ) 

= Et/?(m) 


( 2 . 2 . 8 ) 


where, 


1 \ 


ip{rn) = 


(2.2.9) 




E, the matrix of system coffecients, is the Sylvester matrix associated with the 
polynomials -B{q~ l ) and A(q~ l ). So if A^ -1 ) and B{q~ l ) are coprime, then E is 
nonsingular, and vise versa. 


R = Exl>(m)ift T (m) = ~R3 T (2.2.10) 

and 

R = E4im)^ T {m) (2.2.11) 

So, R is positive definite if and only if E is nonsingular and R is positive definite. 

If the system is unidentifible, which implies A(z) and B(z) are not coprime , then 
we can not find the input which will guarantee to persistently excite the system. 
Before going further, we will use the following definition: 
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Define f(jfc) = ( V T {k - 1), V T {k - 2), , V T (k - L)) T . V{k) is said to be persis- 
tently exciting(p.e.) of order L if there exist Qi, a 2 > 0, such that for all / , there 
exist m, corresponding to which the following equation holds, 

/+m 

( 2 . 2 . 12 ) 

k=t 

It is easy to show that R > 0 if and only if V(k) is persistently exciting of order 
2 n — 1. Since the V(k) are composed of frequency coefficients of the input signal u(f), 
it means that input signal should have a wide frequency band. 

So for noise free systems, R > 0 if and only if A(z) and B(z) are coprime and V(k) 
is p.e. of order 2n — 1. 


Case 2. e(m) ^ 0: 


Assume that the signal and noise are not correlated, i.e. 


Let 


Then 


R = E 


( Z(m-l) \ 

Z(m — n) 
V(m — 1) 

V V / (tt7 - n + 1) / 


EV(i)e T {j) = 0 V ? , j 


m = 

e(k) = 


-A^r ik) 


(2.2.13) 


(2.2.14) 


( Z T (m — 1), . . . , Z T (m — n), V T (m — 1), ... , V T (w — n) ) 
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( e(m - 1) \ 


+E 


e(m — n ) 


0 


( e T (m - 1), . . . ,e r (m - n),0, . . . ,0 ) 


/ Rz Rzv \ /Re 0 
\R T ZV Rv I l o 0 


V o / 


It is obvious that the necessary condition for 72 > 0 is Rv > 0. If we assume 
R e > 0, using the argument of algebra, it is easy to prove that Rv > 0 is the suffcient 
condition for R > 0 as well because Rz ^ 0. 


So assuming R e > 0, then R > 0 if and only if V(A:) is persistently exciting of 
order n — 1. 


The above analysis is very useful for system identification because the condition 
of persistent excitation of V(k ) only depends on the input signal u(fc), not involving 
any output signal y(t). For the nonlinear polynomial systems , quite same argements 
can be advanced and similar results will be obtained. We omit the analysis here. 


2.3 A Least Squares Algorithm for System 
Structure Estimation 


In this section, we consider to use the least squares algorithms to identify system 
parameters and estimate the order of a system. The idea is to minimize the mean 
square error by adjusting the system parameters and order based on the prediction 
error modeling as depicted in Fig. 2. 1(a). 
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Suppose we choose n as the order of the system. With reference to equ.(2.2.3), we 


have: 


W(m) = -a„ilf(m - 1) — a n2 W(m - 2) - ... — a nn \V(m — n) 
+b n0 V(m - 1) + . . . + V(m - n) + e n (m) 


N > m > n 


(2.3.1) 


Define: 


Then, equ. (2.3.1) becomes: 


X(k) = 


W(ifc) 


(2.3.2) 


[W r (n),...,^(AT)] = [^ nil ,...,^ n ] 


X(n-l) X(n-2)... X{N- 1) 
X{n- 2) X(n-l)... X{N- 2) 


X(0 ) X(l)... X(N — n) 


+ [e„(n),e n (n + 1), . . . , e n (JV)j = 9^1 + el 


where 


@n,j — ( ^nj- 1) 

0 T n = [0n,lA,n-2,... A,n] 

= [e„(n), e n (n + 1), . . . , e n (A r )] 
(X{n- 1) ... X(N- 1) > 


r T = 


X(0) ... X(N-n). 


(2.3.4) 


The parameter vector 9 n is determined by minimizing a weighted least squares 


criterion: 


e„ = Y. 
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By some algebraic calculation, the LS estimation of 9 n can be specified by : 

k = [rJrJ-'rJW.tAT) (2.3.6) 

W n (N) = [W», W(n + 1), . . . , H''(A')] T (2.3.7) 

To determine an appropriate order of the model, the parsimony principle will be 
a very useful rule. It says that given the input and output data, if several dynamic 
models fit the data well, then the simplest model should be chosen, i.e., the model 
with the smallest number of independent parameters will be desired. Consider the 
situation here, we have a sequence of model structures of increasing dimensions. Given 
observed input and output data {u(t), y{t)}, a better fit will be obtained if we increase 
the complexity of the model structure, i.e., increasing the order n. The essential thing 
is to investigate whether or not the improvement will be significant in some sense . 

The problem considered in this context is using E n as the criteria. The fact is 
quite obvious. If the true order of the system is n 0 , then for n < n 0 , E n will decrease 
significantly for n increasing. When n approaches 7io . the rate of decrease of E n w T ill 
be slow and E n will almost remain the same value for n > n 0 ■ The typical function 
E n of n looks like Fig. 2. 1(b). 

Assuming a maximal system order is No ■ we choose n as an estimation of n based 
on: 

h — min{n | D n < 8 , 1 < n < Ao} 


where i 0 could be choosen as 1 ~ 3, and 8 > 0 is a preselected threshold. 

Also for on-line identification and for decreasing computational burden, the lat- 


En — E n +{ 

i=l 


E , 


(2.3.8) 
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tice filtering structure could be employed because of the Toeplitz structure of the 
coefficient matrix in equation (2.3.3). See [Sode89] for more detail. 
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2.4 Computer Simulation Results 


2.4.1 Some Computational Considerations 

There are several quantities associated with our identification scheme that need to be 
chosen before setting up the algorithm. These include T , M and N. [0, T] is the finite 
time interval of observed data. M is associated with maximal frequency components 
of the input and output signals. N is the maximum order of the sinusoid modulating 
functions. 

For linear systems, we understand how the choice of T and M will affect the 
identification procedure, thus, u> 0 = 2n/T essentially plays the role of a resolution 
frequency. See [PearSSb]. From a heuristic notion, u,’ 0 should be small enough to 
distinguish the characteristic modes of the system. Also, if the resolution frequency 
is small enough, we can choose M large without involving much high frequency 
measurement noise. Then the available data in our discrete autoregressive model 
(2.3.1) and (2.3.3) will increase and the estimation (2.3.6) will become much more 
robust. All this means that T should be large, i.e., we should have enough available 
input and output data. For linear system identification, the choice of T and M can 
be guided by the frequency domain consideration, i.e., the frequencies retained in 
the input and output signal (U, T) should cover the system bandwidth. The highest 
frequency coefficient in the modulating function Mu 0 should be comparable to the 
system bandwidth W c . Also, the Fourier based modulating functions act as a filter. 
Too large Mu q sometimes is not desirable because it might give undue amplification 
to the high frequency noise present in the input and output signals. 
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Assuming that we know the bandwidth of the system to be approximately W c , 
a reasonable selection of ujq and M could be Mlo o ~ 1.25VF C , or a litter bit larger. 
For nonlinear systems, the characteristic modes and system bandwidth are not well 
defined concepts. But there do exist some intuitive physical meanings which relate 
to above discussion. 

The choice of N is very flexible . Here we let N = 2 M to generate a scalar 
discrete autoregressive model. That will simplify the calculation of the identification 
algorithm. 

The computational cost before constructing the discrete model for parameter es- 
timation will be the integrals of the input and output data. Using the Fourier based 
modulating functions not only gives us the direct frequency domain interpretation of 
the algorithm, but also provides us with an efficient computational tool via the FFT 
technique. 

Let z{t ) denote an observed input or output signal on the time interval [0,T’]. The 
basic computational costs are the following integrals for z(t): 

F z{t)e jmuJot dt , m = 0, 1, . . . ,M 

Jo 

First we can sample z(t) uniformly by sampling interval h in generating the dis- 
crete signal Zi = z(ih), h = T/Nj , i = 0,1,..., Nj . Then the integration can be 
approximated by the standard parabolic rule : 

f z(t)e jm “ Dt dt = |[ z 0 + z Nj + 4 £ z ' Wm ' + 2 E 2 ‘ M/ ’ rn ' +°( h4 ) 

Jo t=2,4,... 

m — 0, 1 , . . . , M 
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where W = e N > and o(h 4 ) is the order of the error as a function of the sampling 
interval h . 

Choosing Nj as a power of 2, neglecting the high order error term, the right-hand 
side of above equation can be evaluated by using the usual FFT algorithm, yielding 
the Fourier series coefficients for m — 0, 1, 2, ... , N j-\ : 

Z = \fFT[zv + z N/ , 4zi, 2 z 2 , . . . , 4^-i] 

The computational savings of this algorithm for large Nj are significant , approx- 
imately log 2 Nj/Nj . However, considering the fact that only M Fourier coefficients 
are needed and Nf should be chosen much larger than M for good accuracy in the 
approximation, we need a special FFT-type pruning algorithm. The efficiency of such 
a devised partial FFT algorithm can be demonstrated to be log 2 M/M. 


2.4.2 Some Numerical Examples 

Example 1 . First we consider a low pass system defined by the Chevbyshev filter 
with bandwidth approximately 0.8[rad/sec.]. The transfer function of the system is: 

H( x ^0438 

s 4 + 0.6192s 3 + 0.6140s 2 + 0.2038s + 0.0492 

The objective is to identify the five unknown system parameters 0 = [0.0438, 0.6192, 

0.6140, 0.2038, 0.0492] and to estimate the order of the system. 

Fig.2.2(a,b) show the input/output data on a 40 sec time interval with output 
signal contaminated by 10% RMS additive white noise. The Fourier series coefficients 
for the first M modes are calculated using the first M components of the parabolic 
approximation with Nf — 1024. 
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For comparing the estimated parameters to the true values and checking the per- 


formance of the algorithm, we define the following normalized error criterion: 

I a* 1 = [^E(^) 2 ]>ioo% 

where K is the number of parameters of the system and 0 t , 6 t are the true and 
estimated parameters. Choosing M = 15 and the threshold S = 0.1, a plot of D n 
verses n, and plots comparing the frequency response of the original system and 
estimated system are presented in Fig. 2. 3 and Fig.2.4. The identified system is shown 
below: 

0.0451 

s 4 + 0.6122s 3 + 0.6210s 2 + 0.2167s + 0.0501 

In the low SNR situations, the algorithm will not work as good as in the previous 
case. But if we can get more information about input and output signal, i.e., choosing 
a larger measuring time interval T, the algorithm will give good results as well. For 
the Chevbyshev system, Table 2.1 shows the performance measured by | A# j with the 
different white noise percentage involved in the output signal and two different time 
intervals. From the table we can see that increasing the time interval can improve 
the accuracy of the estimates. 



TABLE 2.1 

LEAST SQUARES ESTIMATION OF THE CHEVBYSHEV SYSTEM 
White noise percentage Measuring time interval T Normalized error | Ad 


0 

[0, 40] 

1.245% 

5 

[0, 80] 

1.835% 

10 

[0, 40] 

3.245% 

20 

[0, 40] 

9.378% 

20 

[0, 80] 

4.113% 
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Example 2. 


Consider a bandpass system with the transfer function: 

_ s 5 4 0.994s 4 + 0.650s 3 + 0.210s 2 4 0.0446s 4 0.004 

~ s 7 4 3.238s 6 4 5.0456s 5 4 4.9335s 4 4 3.1943s 3 4 1.3146s 2 4 0.3063s 4 0.0304 

We use the same binary input signal as used in example one to excite the system, the 
output signal with 10% additive white noise is shown in Fig. 2. 4. Choosing M = 20, 
N j — 1024 and N = 2 M. The quantity D n as a function of n is illustrated in Fig. 2. 5. 
In this case, if we select 8 = 5%, we can obtain the correct estimate h = 7, the exact 
order of the system. Fig. 2. 6 gives the comparison of the magnitude response of the 
original system and estimated system. The identified system is presented below: 

„ _ 1.056s 5 4 0.921s 4 4 0.621s 3 4 0.243.S 2 4 0.0452s 4 0.0043 

~ s 7 4 3.4501s 6 4 5.0903s 5 4 5.100s 4 4 3.067s 3 4 1.3406s 2 4 0.3281s 4 0.031 

For the band pass system, Table 2.2 presents the normalized estimation error with 
various noise to signal percentages and two different time intervals. 


TABLE 2.2 

LEAST SQUARES ESTIMATION OF THE BAND PASS SYSTEM 


White noise percentage 

Measuring time interval T 

Normalized error | A 6 

0 

[0, 40] 

2.567% 

5 

[0, 80] 

3.223% 

10 

[0, 40] 

6.297% 

20 

[0, 40] 

19.473% 

20 

[0, 80] 

11.468%, 


As in example one, increasing the time interval can improve the estimation accu- 
racy and reduce noise effects. Also, based on several simulations, we have found that 
the best selection of 6 is around 0.05 ~ 0.1 which gives quite good order estimation. 
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t (sec.) 

Fig. 2. 5(a) Output signal y(t) contaminated by 10% white noise. 



n 


Fig. 2. 5(b) D n verses n. 
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Chapter 3 


Frequency Analysis Using Short Time 
Transient Data 


3.1 Introduction 

In this chapter, we present some algorithms for frequency response estimation by using 
short time transient input-output data. The development of the complex sinusoidal 
modulating functions from section 3.2 to section 3.5 follows [PearAE, Pear91], and 
some of A. E. Pearson’s unpublished notes. 

In addition to direct sinusoidal steady state measurements, available algorithms 
for frequency response determination are the direct DFT ratio, and the classical cross 
correlation method[Astr75, Ljun87, Sode89]. If a system is in the rest status initially, 
the length of available input and output data is long enough , and the tail part of the 
system impulse response vanishes dramatically, then the frequency response of the 
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process can be estimated by the discrete Fourier analysis: 


G(M) 


Y_ (M) L 
U(MV 


1,2,. ...M 


where, 


_ _ / .27 Th ^ ^ _ 2ttiAt At 

l 7 0-7“ ) ~ At E y(*A<)e J i 

^ i=0 

TT ..2ltk' 0^ , . x 

U(j-j-) & At 2J u(iAt)e J L 
L t=0 

it = 1,2, . . . , M 


(3.1.1) 


(3,1,2) 


At is a sampling interval, Z, is the length of input-output, data and has to be large 
enough to make equation(3.1.2) approximately hold. 


If input-output data satisfies the stationarity assumption, we can use the cross 
correlation method to estimate the transfer function by the following formulations. 


G{juk) — 


RyuiM) , , n 

7WM)’ 1,2 


,.,M 


(3.1.3) 


where, 


and 


2nk 



E r yu (r)e^ 

L/ 

T=~Z/+1 

^( i2 : k ) = 

E r Utt (r)e^ 

L 

+ 

I 

II 


^ L— 1— max(r,0) 

r y u{r)=- E y{At(i + T))u{Ati), r = 0, ±1, ±2, . . . 

«=— min(r,0) 

| L — 1— max(r,0) 

r uu (r) = — E u(At(i + t)u(AU), r = 0,1,2,... 

L t=0 

r* x ,(-r) = r uti (r ) 


(3.1.4) 


(3.1.5) 
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Using these approaches to estimate the spectral densities sometimes will give a 
poor result even if noise is neglibible. The problem is that the error in calculating the 
spectral estimates does not approach zero (in the mean square sense) as the number 
of data points tends to infinity. 

A critical situation is that the model is not at rest initially and the available length 
of input and output signals is not large. Also the signals involve sensor noise. In this 
case, the approximations will entail an even larger error. 

Here, we deal with the general situation when the system is not at rest initially 
and only short transient data is available. The data consists of the input and output 
signals over a set of finite time intervals {[f,, t, + 7 1 ], ? '= 1,2,. ..,£?}• The length of 
each time interval is T = and it is not necessary that these intervals be disjoint. 

[Pear91] proposed a method which used a least squares algorithm to estimate the 
frequency response at selected frequencies by only using short time transient data. 
The advantage of the approach is that there are no assumptions of statistical station- 
arity of the data or steady state operation. In order to get good estimates, the input 
and output data should contain enough energy at the selected frequencies, or a reg- 
ularized algorithm might be employed. The basis of the technique is the modulating 
functional approach. We have constructed the complex sinusoidal modulating signals 
for parametric frequency response estimation. Several formulations of the standard 
normal regression form for the least squares estimation of frequency response related 
parameters have been developed. 

Here a review of the MFFM(Modulating Function Frequency Method) is presented. 
An algorithm of applying the MFFM to MIMO systems frequency analysis is also 
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proposed. Some computer simulation results will demonstrate the performance. 


3.2 Complex Modulating Signal Construction 
and the Modulating Property 

We consider a stable linear system with a single input u(t) and single output y(t). The 
system will be expressed in the differential operator format involving some modeling 
errors and noise: 

A (p)y(t) = B(p)u(t) + e(t) (3.2.1) 

where, 

A(p) = a n p n + a n _ip n 1 flip 1 + Go 

B(p) - b m p m + b m -ip m 1 + • • • + b\p l + b 0 (3.2.2) 

A(p) and B(p) are polynomials in the differential operator p = d/dt , m < n. Modeling 
errors and measurement noises will constitute the term e(<). Given the input-output 
data [u(t), y(t )] over the finite time intervals { f.-t-T], i = 1,2, . . . , Q }, the problem 

is to estimate the actual frequency response G(ju>) = B{ju>)j A(ju>) at the frequency 
knots {A’cjo, k = 0, 1, 2, . . . , A/}, where ujo = 2 it /T is a user selected “resolving” 
frequency and M a user chosen positive integer. The time intervals -f T ] need 

not be disjoint, but experience suggests that twenty to fifty percentage independence 
will be necessary in order to avoid singularity of the coefficient matrix for the normal 
regression equations in the LS estimation. Generally, the degeneracies in the LS 
estimation will not occur in normal operating records if the time intervals are more 
than half disjoint. 

In chapter 2, we have discussed that 4>{t) is a modulating function of order n over 
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a fixed time interval [0,T] if 4>(t) is sufficiently smooth and satisfies the end point 
conditions: 

^•)( 0 ) = ^(0 (T) = 0, i = 0, 1,2, . . . ,n - 1. 

Clearly, many modulating functions could be constructed which satisfy the condi- 
tion. Here a specific set of complex sinusoidal modulating signals are built which is 
conducive to solving the problem at hand. 

Define: 

<t> m {t) = e jmwot {l - e ]UJot ) n , 0 < t < To = ^ 

m = 0, 1, 2, . . . , M (3.2.3) 

Using the binomial expansion (a + b) n = ELo C k a k b n ~ k and changing the index 
of summation, we can equivalently rewrite (3.2.3) as 

<Pm(t) = £ b k e* m + k) ^ (3.2.4) 

Ar=0 

where b k relates to the binomial coeflBcient: 

h = (-1 ) k c k 

Obviously, equations(3.2.3-4) define a set of modulating functions of order n over 
a time interval [0,T]. One of the important motivating factors for building these 
complex modulating signals is that modulating the input and output data by these 
signals will entail calculating the Fourier series coefficients at the frequency knots 
km o, k = 0, 1, 2, . . . , M -f n and automatically build some algebraic equations involving 
the estimation of the system transfer function. This is the important modulating 
property and will be discussed further below. Another factor is that equations(3.2.3- 
4) comprise linear combinations of commensurable complex sinusoids and the FFT 
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algorithm can be applied to compute the Fourier series coefficients of input-output 
data over the time intervals [t,,t, + 7]: 


Z k (i)= f z(t + t t )e~ M di 
Jo 

= 0, 1, 2, . . . , M + n, j = 1,2,...,<2 


(3.2.5) 


or in terms of sine and cosine transforms: 

r T 

Ztii) = / z(t + t{) cos ku> 0 tdt 
Jo 

Zt(i) = f z(t + t,) sin kuJotdt 
Jo 

These Fourier transform series coefficients can be computed accurately and efficiently 
because of the availability of the parabolic approximation and the FFT technique(see 
chapter 2). 


Modulation Property 


Let P(p) be a differential operator of order at most n, i.e., a polynomial in p — d/dt 
of degree < n, and z(t) any sufficiently smooth function defined on [0,7’]. Then the 
modulation of P(p)z(t) with <j> m (t) over [0,7] will satisfy 

»T m+n 

/ <t> m (t)P{p)z(t)dt = h- m P{-jkuJ 0 )Z k (3.2.6) 

J ° k—m 

where Z k is the k th harmonic Fourier series coefficient of z(t), i.e. , 

Z k = [ T z(t)e ]kuot dt 

Jo 


Proof: To prove the modulation property (3.2.6), we only need to verify the special 
case of P(p) = p' for a fixed integer i < n. Due to the superposition property of the 
polynomial, the result for a general n th degree polynomial will follow immediately. 
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For sufficiently smooth function z(t) defined on [0.7 1 ], the left side of (3.2.6) in this 
special case is 

f <f> m (t)p' z{t)dt = (-1)' ( z(t)p l 4> m (t)dt (3.2.7) 

Jo Jo 

Equation(3.2.7) is obtained by integration by parts n times over [0, T] taking into 
account the end point value properties of the modulating functions. The crucial point 
is that none of the initial point derivatives z*’*(0) and z^(T) will appear because of 
the end point conditions. Substituting the expression (3.2.4) of into (3.2.6), 

using equation (3.2.7) and changing the index of summation, verifies the modulating 
property for P(p ) = p ' . 

In the following, we formulate several least squares algorithms to estimate the 
system frequency responses at the frequency knots kuj 0 , k = 0, 1 , 2 , . . . , M + n using 
the complex sinusoidal modulating functions. 


3.3 Least Squares Estimation ( The First 
Formulation) 


Rewriting the linear model (3.2.1) in the equation error form, modulating both sides 
by the modulating signals <?!> m (f), we can obtain: 

[ (j> m (t)[A(p)y{t) - B(p)u{t)\dt = [ <f> m {t)e{t)dt (3.3.1) 

Jo Jo 

Based upon the modulation property (3.2.6), the preceding equation could be repre- 
sented by the form: 

m + n m+n 

£ b k - m [A(-jku 0 )Y k - B(-jkuo)U k ] = £ h. m E k (3.3.2) 

k=m k-m 
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Define the real and imaginary parts of the polynomials A(jku>o), B(jku)o),Y k .U k 

by: 

A{jkw o) = a k + j 0k 

B(jku> o) = 7jt + j&k (3.3.3) 

Also assume: 

n = n e +;n - 

U k = U' k +jU‘ k (3.3.4) 

Rewriting equation(3.3.2) in terms of the real and imaginary parts, we can obtain 
the following equations: 

m + n 

£ + An* - 7*e* c - Ai'ii = s” 

k~m 

m+n 

E [-o*n* + An* + lit? - «£] = 4 (3.3.5) 

k=m 

where, 

m+n 

4 = *(£**—£'*) 

k—m 
m + n 

e I m = %Y. b ^-mEk) 

k=m 

Computing the Fourier series coefficients Y k (i) and U k (i) corresponding to each 
time interval [f,,/, + T], i = 1, 2, . . . , Q, k = 0, 1, . . . , M by the equations: 

Y k (i) = Fyit + W^dt 

J 0 

n(0 = y t c (i) + jY((i) 

Uk(i) = f u(t + u)e’ ku ° , dt 
Jo 

Uk(i) = U‘ t (i)+ jV t (i) (3.3.6) 
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tM>) = 


and define a 2 x 4 data matrix V’r(r) in terms of these harmonic series by : 

■ mo mo -mo -v‘ t (iy 

-Y{(i) Yifi) U‘ k (i ) -tfj(i). 

Define a 4 x 1 parameter vector by : 

' Oik' 


(3.3.7 


Ok 


0k 

7 k 

Ok J 


(3.3.8) 


Then equation (3.3.5) could be represented by the 2 x 1 vector equation form as 
follows: 


m+n 


r ^(*) 

LdwJ 


Y bk- m tl’k{i)0k = 


A:=m 


— £m (0 


m 


= 0, 1, 2, . . . , M, i = 1,2,.. .,Q 


(3.3.9) 


For each fixed i, we can collect all the 2 x 1 vector equations for m = 0, 1, 2, ... , M 
and form the following regression equations for serving the least squares estimation: 


♦(*)*-£(*)=£(*') 


i = 1,2,- ..,<3 


(3.3.10) 


where. 


*(») = 


Uo(i) b 2 ip2(i) ••• b n ipn(i) 0 

0 bo’l/ , i(i) b\ Ip 2 ( i ^n— 10n+l(O b n tpn+i^i') 


0 

0 


L o o o 


••• b n yj n+M {i ). I 

(3.3.11) 
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6 = 

‘“To 
0 1 
0 2 

» £(*) = 

'»/(*)' 

0 


- - 


.0 . 


i( 0 = W) 

Here we assume that ao = 1 without loss of generality. 

Collecting these equations for i = 1,2, forming the standard regression 

equation form and choosing the 0 to minimize the least squares criterion: 

min e = min(<l>0 — () T (&0 — £) (3.3.12) 

6 6 

where, 



'«*>(!) ' 


■f(l) ' 

<t> = 

*(2) 

II 

c 

<3 

e(2> 


.*(Af). 


MM). 


There are no disjoint requirements to these finite time intervals {[f,, f, + T],i = 
1, 2, . . . , Q}. Of course a certain degree of independence in the input and output data 
over these time intervals is needed. The input and output data must contain a suffi- 
cient energy at these frequency knots in order to avoid the degeneracy of solving the 
least squares estimation (3.3.12), i.e. , the matrix $ being composed of the harmonic 
Fourier series coefficients of the input and output data should have full rank. 


Also, one kind of “ bootstrapping" least squares algorithm can be developed by 
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stages according to the values assigned to m in order to avoid degenarices of I 1 LS 
algorithm due to the large number of unknown parameters. 

Initial Stage (m = 0): 


Yo(0 = b kipk(iWk 

k=0 

?= 1,2,3,..., A r 


(3.3.13) 


where the data-related quantities are defined by the following vector-matrix equations: 

f-Ya C (i)\ (-m)\ 

>b(0= j V*(0=l o J ( 3 -314) 

The definitions of ^(i) and 9 ^ for k > 1 are exactly the same as (3.3.7) and (3.3.8). 
Oo = 70 - 


Stage m (m = 1, 2, . . . , M): 


Again a set of 2 dimensional vector equations are derived using (3.3.9) for the case 

of 777 > 1: 

Y m (l) = M„+m(i)0n+m + €m(0 (3.3.15) 

i = 1,2 ,. m = 1, 2, . . . , M 


where Y m (i) is defined in terms of the parameters and data of the preceding stages 


by: 

n+m-1 

Y m (i) = - £ b k - m j>k(i)0k 

k~m 

The residuals for all stages are given by: 


n + m 

~ ^ bk—m 

k~m 



(3.3.16) 
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(3.3.17) 


From the equations(3.3.3) and (3.3.8), it is easy to see that the transfer funtion 


G(jku o) = B(jku)o) / A(jkujo) relating to the parameter vector 6k is given by: 


X(G(jkul o)) = + JT- 

$(G(jkuo)) = 

a k + Pk 

or equivalently by the magnitude- phase relation: 


\G(jku> 0 )\ 2 = 


<*1 + ft 


Z.G{jktjJo) = tan 1 — — tan 1 — 
7fc a *r 


(3.3.18) 

(3.3.19) 


(3.3.20) 

(3.3.21) 


Investigating this least squares estimation procedure leads to the following remarks: 


1. The frequency knots k.us 0 , k = 0. 1, 2, . . . , M-f n scatter in the frequency range 
from 0 to ( M + n)ujQ. Hence, if we know the bandwidth W c of the system, 
then a reasonable selection of uio and M should satisfy Mu>q ~ 1.25iT c , or 
a little bit larger if it is desired that the estimation of the system frequency 
response cover a range 25% larger than the bandwidth W c . 


2. For bandpass or high pass systems, <j> m (t) could be modified by: 


= e j{m+mo )““'(! 


- e JU,ot 


) 


m = 0 , 1 , 2 , . . . , M 


In this case, the transfer function estimate covers the frequency range from 
uiqloo to (M + n T mo)wo- We are only interested in the frequency band above 
the frequency moUo . The available information about the system cutoff 
frequency and bandwidth could be utilized to determine the choice of m 0 and 
M. 
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3. The total number of unknown parameters, i.e., the size of 9 is 4(M + n). 

Since the dimension of each 4> is 2Q(M + 1) x 4(n + M), the total number 

of time intervals + T ] should satisfy Q > . For good accuracy, 

the number of equations should be double the total unknown parameters, i.e., 

O > 4 ( M+W ) 

hf - (M+l) ' 

4. If the resolving frequency u?o has to be small for some special practical appli- 

cations, M should be large in order to cover the system bandwidth. In this 
case, the total number of unknown parameters 9k is incremented. The matrix 
4> r 4> tends to be singular and the least squares solution for the eqn.(3.3.12) 
will be degenerate. There are two methods to solve this problem. One is to 
try obtaining extra available data, i.e., increasing the total number of time 
intervals. Sometimes this is expensive. Another one is to exploit a numer- 
ical advantage of solving eqn.(3.3.12) which is based upon the fact that the 
equations(3.3.9) and (3.3.10) are partially decoupled with respect to the in- 
dex m. Parameter 9 k only depends on the 9 t for i < k. Therefore, one kind 
of “bootstrapping” of the least squares algorithm could be utilized to divide 
estimating the parameter 9 into M+l stages. • • • • 9 n ) can be identified 

at the first stage, corresponding to m = 0, which requires Q satisfying N > 2 n 
for 4 n unknown parameters. For good estimation, Q should be larger than 
4n. For each succeeding stage (m = 1,2,...,M), only 9 n+m (4 unknowns) 
needs to be estimated, which only requires Q > 2. By this way, the number 
of unknowns at each stage is kept to the lowest level and most likely the 
degeneracy problem will be prevented during the least squares procedure. 

5. Noting the fact that |0*| goes as k n as k increases for k = 1, 2, . . . , M + n, a 
nonlinear transformation might be necessary to scale 9 appropriately for good 
numerical accuracy, or we can normalize all parameters by pVF”, where p is 
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a scalar constant and W c is the system bandwidth. 


We have developed a scheme to estimate the original transfer function G(jku) o) 
at frequency knots koj 0 ,k = 1,2, . . . , M + n. The system does not have to be in 
the rest status initially and the input and output data could be time limited and 
transient. The data matrix <pk{i) in equation(3.3.7) comprising the k th harmonic 
frequency series in the i ih time interval, k = 1, 2, . . . , M + n, i = 1,2, ...,<5, is 
mutually orthogonal. This reveals the maximum utilization of information contained 
in the input and output data for our least squares procedure. Clearly it is a direct 
result of the Fourier nature of the underlying formulation. 

One numerical problem in our scheme might arise due to the highly nonlinear rela- 
tions in the equations(3.3.18-21) which involves the difference between the estimated 
parameters whose values may be large for large frequencies. In the next section, 
we will develop other schemes, which have different structure from this one and will 
alleviate the potential numerical error source mentioned above. 


3.4 Least Squares Estimation ( The Second 
Formulation) 


For the original model described by equation(3.2.1), multiplying both sides by A(— p), 
we can see that the input and output data also obey the following relation: 

A{-p)A{p)y{t) = A(—p)B(p)u(t) + A(-p)e{t) (3.4.1) 


42 


Utilizing the complex sinusoidal modulating signals of order 2 n defined by: 


<f> m {t) = e jmLJ0t (\ — e 3ulo1 ) 2n 

m + 2n 

£ 

k—m 


k—m 

to project the model into the algebraic equations, we get the result: 


(3.4.2) 


m+2n 


where, 


h-mA{jku} Q )[A{ -j ku 0 ) Y k - B{-jkio 0 )U k ] = e m 

k—m 


h = {-i) k cl 

m-f 2n 

£m = h-mA{jkuo)Ek 

k—m 


(3.4.3) 


Define: 


A(jkuJo)A(-jkiJo ) - ak 
A(jkuJo)B{-jkLOo) = a* + j0k 


(3.4.4) 


(afc,Qfc,/3jt) are real numbers. 


Rewriting (3.4.3) in terms of real and imaginary parts of the equations, we obtain: 

m+2n 

£ 

k=m 


£ h. m \a„Y£ - a t Ui - Mi] = 4 

m+2n 

£ 5*-4 a„Y; - a t Ul - Mi] = 4 

k — m 


(3.4.5) 


Define a 2 x 3 data matrix fl'fc(i) involving the k th harmonic series of input output 
data in the i th time interval [<,, f, + T\ by : 

■n e (0 -w) utw ' 

.n s (0 -urn) -uzd). 


®fc(0 = 


(3.4.6) 
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Define a 3 x 1 parameter vector by: 


Qk 


6 k — 




LAJ 


(3.4.7) 


Then we combine the real and imaginary parts of equation(3.4.5) into a 2 x 1 
vector equation as follows: 

m + 2n 

^ 1 ^k—m^k^k — £m(^) 
k=m 


m = 0, 1 , . . . , M, i = 1 , 2, . . . , Q 


(3.4.8) 


Collecting all these equation and rewriting it into a standard regression format, 
which is analogous to eqn.(3.3.12), the least squares algorithm can be employed to 
estimate the parameters (0 U 0 2 , . ■ . , 0 M+ 2 n)- In this case, the relation between the 
estimates of the parameters and the real and imaginary parts of the original transfer 
function are: 

*(GO'*wo)) = — 
ak 

Z(G(jku>o)) = (3.4.9) 

&k 

or equivalently by the magnitude and phase relations: 

\/oti + 

\G(jku, 0 )\ = ^ — — 

lG(jku 0 ) = -tan _1 — (3.4.10) 

a* 

We make the following remarks by considering the implementation of this scheme: 


1. The estimates of the frequency response of the original system cover the fre- 
quency band from 0 to ( M + 2n)u; 0 . An approprite choice of the parameter 
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M and ujq should satisfy the relation: 


(M + 2r?> 0 ~ 1.25 W c 

2. The total number of unknown parameters in this case is 3(A/ + 2n). Hence, the 
number of time intervals Q should satisfy 2 Q(M + 1) > 3 (M +2n). However, 
in order to avoid the degeneracy of the least squares estimation, the number 
of equations should be double the unknown parameters, i.e., Q > 

3. Employing the partially decouped nature of equation (3.4.8) might be nec- 
essary to solve the least squares solution if a large number of unknown pa- 
rameters need to be estimated. In this case, the number 6n of unknown 
parameters (0j, 0 2 , ■ ■ ■ , ^2n) are involved in the first stage corresponding to 
m = 0. The choice of Q should satisfy Q > 6 n. For each succeeding stage 

(m = 1,2, M), only 3 unknowns need to be estimated assuming the use 

of the preceding estimate of the parameters. Comparing to the first formula- 
tion presented before, the total number of unknowns at the first stage is 6n, 
verses 4n unknowns. But only 3 unknowns are required to be estimated in 
the succeeding stages, verses 4 unknowns for the first formulation. 

One point we should mention is that the equat,ions(3.4.9-10) avoids differencing 
the large quantities associated with the estimated parameters to calculate the fre- 
quency responses at high frequencies. This is a numerical advantage of the second 
formulation. 
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3.5 Least Squares Estimation ( The Third 
Formulation ) 

We can develop a dual formulation to the previous one. Multiply the original model 
(3.2.1) by B(—p) shows that input output data also should satisfy the relation: 

B(-p)A{p)y(t) = B(-p)B(p)u(t) + B(-p)e(t) (3.5.1) 


Again using the complex sinusoidal modulating functions of order 2 n to project 
model (3.5.1) into the algebraic equation: 

m + 2n 

H h m B{jkuio)[A(-jkuio)Y k - B{-jkuj 0 )U k \ = e m (3.5.2) 

k=m 

m = 1,2, 3, . . . , M 

where, 

m+2n 

£m — ^ ^ ^A: — m J ) Ek 

k=m 

Similarly, define: 

B(jkuj 0 )B(-jku! 0 ) - b k > 0 

A(jku> 0 )B(-jku 0 ) = a k + jfik 

'b k ' 

Ok = a k (3.5.3) 

.—0k. 

Define a 2 x 3 data matrix ^ k (i) for the i th time interval by: 

rV£(i) - Y t ‘(i ) Vi*(0 1 

"MO = (3.5.4) 

Lc£(o -n*(>) -nioJ 
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Then equation (3.5.2) can be rewritten into a 2 x 1 vector equation which is 
analogous to eqn.(3.4.8): 


m 


m+2n 

53 h-m*k(i)9k = em ( 0 
k=m 

= 0,1 i = 


(3.5.5) 


A standard regression from could be constructed for setting up the least squares 
estimation for the unknown parameters ($i, 9 2 , ■ ■ . , 0m+ 2 n)- 


Calculating the frequency responses in terms of the estimated parameters reveals 
the following relation for the real and imaginary parts: 

c*kbk 


%(G(jku 0 )) = 




(G(jkuo)) — 


<* l + PI 

Pkh 


a ; l + PI 

or, equivalently for the magnitude and phase relation: 

h 


|G(j^„)| = 




i 0k 

/.G{jku> 0 ) = — tan 1 — 


(3.5.6) 


(3.5.7) 


The data matrices defined in equations (3.4.6, 3.5.4) have the same structures by 
interchanging the direction of input and output flow. This reveals the dual property 
between these two formulations. 


3.6 Computer Simulation Results 

Here some typical numerical results are presented to demonstrate the performance of 
the algorithm. We employ the first formulation for the following two simulations. 
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Example 1. First we consider a low pass system with a bandwidth approximately 
12[rad/sec.]. The transfer function of the system is 

H (s) = 1.7s 2 + 1736.8 

~~ s 3 + 19.1s 2 + 257.48s + 1700.0 

A random binary process was utilized as an input signal to excite the system because 
it contains sufficient energy in a broad frequency band. Data was collected in 8 
seconds. The plots of input signal and contaminated output signal with SNR=15(dB) 
are shown in Fig.3.1(a) and Fig. 3. 1(b). The whole data set was divided into the six 
time intervals [t^ti + 3] for ! = 0, 1, ... ,5. The sampling frequency f s was 200(Hz). 
A 1024-FFT was utilized in calculating the Fourier series coefficients on each time 
interval. We estimated the frequency response Hi (jus) at frequency knots {ku>o, k = 
0, 1, 2, . . . , 18}, u;o = 2.09. The frequency knots cover the passband, transition band 
and a part of the stop band of the system. 

Fig.3.2(a-b) shows the results of the frequency response estimation for a noise- 
free case, along with the results obtained by applying the direct ratio of the Fourier 
transforms and the standard cross correlation approaches to the same data. 
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Fig. 3. 2(a) Estimated magnitude response at the frequency knots {fcu> 0) k = 
0,1,2,..., 18} in the noise-free case. 
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Fig. 3. 2(b) Estimated phase response at selected frequency knots in the noise free 
case. 

Also the simulation results are shown in Fig.3.3(a) and Fig.3.3(b) for the ouput 
signal contaminated by white gaussian noise with SNR=15(dB). 



Fig. 3. 3(a) Estimated magnitude response at selected frequency knots for 
SNR=15(dB) case. 


51 



0 


I - 100 


ja -200 
Pl, 




—Trot Phase Response 
ne Our Parametric method 
-eCroa Correlation 
--Direct Ratio rf Fourier Tram 


0 5 10 152025303540 


Frequency(rad/sec.) 

Fig. 3. 3(b) Estimated phase response at selected frequency knots for 
SNR = 15(dB) case. 


From the graphs we can see that noise will not affect the algorithm at low fre- 
quencies. But it does increase the estimation bias for high frequencies. The cross 
correlation and direct ratio methods produce large oscillations. 


Example 2. The second example we considered was a high pass system with the 
transfer function: 


H 2 (s) = 


s 3 + 22.02s 

s 3 + 22.24s 2 + 247.44s + 1943.23 


The same random binary signal was used to excite the system. Again the data 
was collected in eight seconds. Fig. 3. 4 shows the output signal contaminated by- 
white gaussian noise with SNR=20(dB). As in example 1, the whole eight second 
data was divided into the six time intervals + 3], t, = 0, 1, .... 5. The sampling 

time interval At was 0.005. A 1024-FFT was used to calculate the Fourier series 
coefficients. Again we estimated the frequency response H 2 {ju > ) at frequency knots 
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{ kuo , k = 0,1, 2,..., 18}, wo = 2.094. Fig.3.5(a)(b) present the estimates of the 
magnitude response as well as the phase response for the noise-free case along with 
the results in applying the direct ratio and cross correlation methods to the same 
data. Also the results are shown in Fig.3.6(a)(b) for the SNR=20(dB) case. 



Fig. 3. 4 Output signal contaminated by white gaussian noise 
with SNR=20(dB) and collected in 8 seconds. 
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Fig. 3. 5(a) Estimated magnitude response at selected frequency knots in the 
noise-free case. 



Frequency(rad/sec.) 

Fig. 3. 5(b) Estimated phase response at selected frequency knots in the noise-free 
case. 
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Fig. 3. 6(a) Estimated magnitude response at selected frequency knots for 
SNR=20(dB) case. 



Frequency(rad/sec.) 

Fig. 3. 6(b) Estimated phase response at selected frequency knots for 
SNR=20(dB) case. 

In the noise-free case, our algorithm gives almost exact results. The direct ratio 
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shows better performance than the cross correlation. Again, noise will not affect the 
estimates for low frequencies, but it increases the estimation error for high frequencies. 
The cross correlation and direct ratio gave large biases at high frequency knots. 

From the above numerical simulations, we can see that if only short time transient 
data is available, our parametric sinusoidal modulating functional method works bet- 
ter than the classical cross correlation and direct ratio methods. If the system is not 
at rest initially, our method will present much better performance. Here one thing we 
need to point out is that if the order of modulating signal is larger than the order of 
the actual system, our algorithm works well. But if the order of modulating signal is 
less than the actual order of the system, our algorithm will fail to give good results. 


3.7 Frequency Analysis for MIMO Systems 


Assume a MIMO model is given by the following transfer function: 

y{s) = H(s)u(s) (3.7.1) 

where u is a p x input vector signal, y is a p 2 output vector signal, H = (/?,,;) is a 
P2 x Pi transfer function matrix. 


Each transfer function element h t _i(s) is given by: 

• A ) ~ 4,,(S) 

B ti i(s) and Aij(s) are coprime polynomials in s. 


(3.7.2) 


Given input-output data [u(t),y(t)j over the finite time intervals { [t,-, t, -f T],i = 
1,2, ... ,Q }, the problem is to estimate the actual frequency response H(ju;) = 
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( hi'i(ju > )) at the frequency knots {ku>o, k - 0, 1,2, for each (*',/) th transfer 

function, where ujq = 2tt jT. 


Based upon the transfer function structure (3.7.1), the problem can be decomposed 
into P 2 subproblems by processing each row individually. The i th element of output 
y corresponds to the MISO system: 


/ ui(s) \ 


JZi('S) — \kt, 1(^)5 • • • 1 ^t,pi(^)] 


\ Up, (s) / 



>M 5 ) 


U/(s) 


(3.7.3a) 


Let /l,-(s) be the least common multiple of {^4i,i(-s), . . • , ^4., Pl (-s)} and assume 

BM 


B u (s) = 


AiAs) 


A(s) 


Then, adding some modeling errors and measurement noise, the subsystem (3. /.3a) 
can be rew T ritten as a differential operator equation form: 


pi 


Ai(p)yi{t) = B,Ap)ui(1) + c.-(0 

;=i 


(3.7.36) 


Therefore, from now on we only consider the MISO systems without loss of gen- 
erality. We also drop the subscript i for easy notation. 


Hence, consider a stable linear MISO system: 

A{p)y{t ) = J2 Bi(p)ui(t) + e(t) 
1 = 1 


where, 

A(p) = a n p n + a n _ip n-1 + . . . + ai p 1 + flo 
Bi{p) = km,p mi + 6 ; , m( _ip m,_1 + . . • + bup 1 + 6/, 0 , m, < n 


(3.7.4 a) 


Z = 1,2,3,. 


r 


57 


It is to be acknowledged that the pairs (A(p), Bi(p)), l — 1,2, ...,r will not 
generally be coprime. However, this is not an issue for us because it is the ra- 
tios Bi(ju)/A(ju:) that we seek to determine for knots ku o, k = 0,1,..., M, rather 
than the polynomial value (A(jw), Bi(ju>)). This avoids a difficult issue of sufficient 
parametrization and minimality for state space models which in our context is a 
separate issue. 

Let {<j5> m (f)} t> e the set of n th order modulating functions defined by (3.2.3). Rewrit- 
ing the linear model (3.7.4a) as the equation error form, modulating both sides by 
the signal we can obtain: 

m+n r m + n 

£ b k . m [A(-jku>o)Yk Wu] = E h k-mE k (3.7.46) 

k=m 1=1 k=m 

Define the real and imaginary parts of the polynomials A(jkuo), Bi(jku^),Y k .Ui f k 
as: 

A{jku > 0 ) = ctkA j@ k 
Bi{jku 0 ) = 7 i,k + j&i,k 

n = n c +j'n*= f T yio^dt 

Jo 

Uijk = U{ k + jUi k = [ T uAty^'dt (3.7.5) 

’ ’ Jo 

Rewriting the equation (3.7.4b) in terms of the real and imaginary parts, we can 
obtain the following equations: 

m+n r 

E + AH - + WJ = *5 

k=m 1=1 

m + n r 

E k-4-«*H + ah + E(-w('u - {,,„[/,%] = e ' m (3.7.6) 

k=m /=1 
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where. 


m+n 

fc=m 

m+n 

53 b k-mE k ) 

k=m 


Compute the Fourier series coefficients Y k (i) and {+*(?') corresponding to each time 
interval [t u U + T], i = 1,2, . . . , Q, k = 0, 1, . . . , M, by the equations: 

Y k (i)= [ T y(t + U)e^ ot dt 
Jo 

n(i) = n c (0+jn*(0 
u,,t( 0 = f u, it + 

Jo 

Ui, k {i) = Ul k {i)+jUl k {0 

and define a 2 x (2 + 2r) data matrix in terms of these harmonic series by : 

Y k c (i) Y£{i) - U{ M {i ) -U{ ik {i) ... -U^ k (i) -U:M 

—Y k (i) Y k (i) Ul k (i) -Ul k (z) ... U; tk {i) - U^ k (i ). 

Define a 2(r + 1) parameter vector by : 

0 k = [a k , 0 k ,~h, k ,6 lik , . . . ,~f T ,k,t>T,k] T (3.7.9) 

Then equation (3.7.6) could be represented by a 2 x 1 vector equation form as follows: 

m+n £m(0 1 

£ bk-mMWk = I = £ m(0 

k=m 


i'k{i) 


(3.7.7) 


(3.7.8) 


Um(0J 

m = 0, 1,2, . . . ,M, t = l,2,... ,Q 


(3.7.10) 


In the same way as SISO systems, we can use a bootstrapping algorithm to form 
a sequence of linear regressions for equation (3.7.10). We omit the algebraic detail 
here. After getting the estimates of 9 k , we can employ formulae similar to (3.3.20-21) 
to estimate the frequency response function for each subsystem Bi(juj)/A(ju) at the 
frequency knots ku > 0 . Other formulations discussed in the previous sections could also 
be applied to MIMO systems. 
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3.8 Numerical Simulation Results 


Here a numerical simulation result is presented to show the performance of the algo- 
rithm. We utilize the first formulation. 


We consider a two input, one output system. The first subsystem is a low pass 
second order system with the transfer function given by: 




2s + 160 
s 2 + 20s + 160 


The second subsystem is a high pass system with the transfer function given by: 


H 2 (s) 


s 2 -f 50s + 54 
s 2 -f 40s + 500 


The input-output relation in the form of (3.7.4a) is 


A{p)y(t) = # i ( p ) u i (*) + B 2 {p) + u 2 (t) 


where 


A(p) = ( p 2 + 20p + 160)(p 2 + 40 p + 500) 

£,(p) = (p 2 + 40p + 500)(2p + 160) 

B 2 (p) = (p 2 + 50p + 54) (p 2 + 20p + 160) 

Notice that the (A(p), Bi(p)) and (i4(p), B 2 (p)) are not pairwise coprime. 


A random binary process was used as an input signal uj to excite the first subsys- 
tem. The superimposed sinusoidal signal 

u 2 = sin(5f 3 ) + sin(20i 2 ) + sin(6f) + sin(25< + 0.9434) 

was used as an input signal to excite the second subsystem. Data was collected in 
15 seconds. The whole data set was divided into the twelve time intervals [ti,ti + 4] 
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for t x = 0, 1,2, ... , 11. The sampling frequency was 200(Hz). 1024-FFT was utilized 
for calculating the Fourier series coefficients on each time interval. We estimated the 
frequency response Hi(jui ) (l — 1,2) at frequency knots {kuJo,k — 0, 1, 2, . . . , 15}. 
w 0 = 1.57. The frequency knots cover a part of the passband, transition band and 
stop band of the two subsystems. 

For comparison, we use ui(f) to excite the first subsystem while holding u 2 at zero 
to get the output and then in a separate simulation use u 2 (t) to excite the 

second subsystem to get the output y 2 {t) while holding u\ at zero. Then we apply 
the direct ratio of Fourier transform and the standard cross correlation approaches 
to the separate data [u^t), yi(t)}, [u 2 {t),y 2 (t)] to estimate the frequency response of 
the subsystems. Although this comparison is favorable to the direct ratio and cross 
correlation, the MFFM still gives better results. 

Fig.3.7(a,c) give the magnitude response plots of the two subsystems by using the 
MFFM for the data [ui,u 2 ,y] in the noise free case along with the results obtained 
by applying the direct ratio and the cross correlation approaches to the data [ui,j/i], 
[u 2 ,y 2 \. Fig.3.7(b,d) are the phase response plots. 

Also the simulation results of the magnitude and phase responses are showm 
in Fig.3.8(a-d) for the output signal contaminated by white gaussian noise with 
SNR=15dB. 

As in the SISO system case, the MFFM gives almost exact results in a noise-free 
situation. Again noise did not affect the estimates of the MFFM at low frequencies. 
Comparing to the SISO case, noise causes larger estimation biases at high frequencies 
for the MFFM. The direct ratio and the cross correlation estimates present a large 
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oscillation and have large error, especially at high frequencies. 



Frequency(rad./sec.) 

Fig. 3. 7(a) Estimat ed magnitude response of subsystem one at selected frequency 


knots in the noise-free case. 



Frequency(rad./sec.) 


Fig. 3. 7(b) Estimated phase response of subsystem one at selected frequency 
knots in the noise-free case. 
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Frequency(rad./sec.) 


Fig. 3. 7(c) Estimated magnitude response of subsystem two at selected frequency 
knots in the noise- free case. 
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Fig. 3. 7(d) Estimated phase response of subsystem two at selected frequency 
knots in the noise-free case. 
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Frequency(rad./sec.) 

Fig. 3. 8(a) Estimated magnitude response of subsystem one at selected frequency 
knots for SNR=15(dB). 



Frequency(rad./sec.) 


Fig. 3. 8(b) Estimated phase response of subsystem one at selected frequency 
knots for SNR=15(dB). 
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Frequency(rad./sec.) 


Fig. 3. 8(c) Estimated magnitude response of subsystem two at selected frequency 
knots for SNR=15(dB). 



Frequency(rad./ sec.) 


Fig. 3. 8(d) Estimated phase response of subsystem two at selected frequency 
knots for SNR=15(dB). 
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Chapter 4 


Schemes for Model Reduction and 
Parameter Identification in the 
Frequency Domain 


The approximation of high-order linear models by lower order systems is a very im- 
portant problem involving a wide range of applications, such as signal processmu' 
and filtering , system synthesis, verification, and controller design. Many researches 
have been done in these and related areas [Sham75, LariSd, MullTO. MahlbU]. I he 
main task of model reduction is to find a simple model structure characterizing tin- 
major behavior of an original high-order system and to simplify the hardware design 
and system analysis. The accuracy of the reduced system really depends upon t lie- 
selection of the error criterion. 

During the last two decades, most of the work utilized the classical approximation 
theory to match the state space realization, Markov parameters ( a part ol the im- 
pulse response sequence ), and etc. These algorithms are based on the classical pade 
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approximation and continued fraction expansion which expand the original system 
into a Taylor series about the origin. The disadvantages of these methods are a low 
accuracy in the high frequency band and the possibility of losing stability for the 
reduced model. 

In order to avoid these problems, some new approaches have been investigated. 
Current methods for model simplification and parameter estimation fall into two cat- 
egories: deterministic and stochastic approaches, each approach adopts either a time 
domain or frequency domain format. Some proposed algorithms in [Mull76, Inou83, 
Wagi86] not only match a finite portion of the Markov parameters but also fit a finite 
portion of the output covariance sequence. For dealing with non-minimum phase sys- 
tem, [Tugn86] presented an optimization algorithm for matching the autocorrelation 
sequence as well as the higher order statistics. Based up the principal component 
analysis and singular value decomposition technique, [MoorSl] gives an algorithm 
to transform an original system into a balanced state space realization and reduce 
a high-order model into a lower one. [HuXi87] combines the Pade approximation 
method and frequency fitting approach to acquire better accuracy in the middle and 
high frequency ranges. 

Most of the investigations mentioned above require some assumptions about the 
original system, or availability of data before starting the algorithm. By deterministic 
means, the usual assumption is that the original model is known in advance. If this 
is not the case, then the first step is to identify it using the measurements of input- 
output signals. By stochastic means, the usual assumption is statistical stationarity, 
i.e., the system is working on a steady status. A long data length is required in order 
to get good estimates of statistical quantities of the original system for acheiving good 
accuracy of a reduced model, especially when the tail part of the impulse response of 
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the original model does not vanish significantly. 


One of main advantages of using the modulating function technique is that we 
can attack continuous time models directly, which is much desired in classical and 
adaptive model reference control. By this means, it avoids the potentially significant 
errors in approximating derivatives from noisy input and output signals, selecting 
a sampling frequency, constructing a transformation for mapping a discrete time 
system to a continuous time system, and etc. All these processes will bring distortion 
to the model. Some are even much more noise sensitive. Another advantage of 
using the modulating function approach is that there is no necessity to deal with the 
complicated initial value problem. 

This chapter focuses on reducing the linear continuous model using only transient 
input and output data and concerns with the task of combining the model identifica- 
tion and simplification process together. Apart from linearity, there are no assump- 
tions of knowing the original model or presuming statistical stationarity of the input 
and output data which are crucial to the previous methods. The approach could be 
realized by two steps. The first step is to set up a least squares algorithm for estimat- 
ing the frequency response of the original model at selected frequency points using 
short input output data records. This issue has been discussed in chapter 3. The 
second step, which we will develop here, is to select a lower order system to match the 
frequency responses at these dominant frequencies. These frequencies should cover 
the system bandwidth in order to characterize the behavior of the model. The orig- 
inal model could be minimum or nonminimum phase and our method can deal with 
phase and magnitude information flexibily. If the phase response is not important, 
we can exploit, this advantage for further reducing the structure of the model utilizing 
an error norm criterion only depending on the magnitude. 
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4.1 Scheme 1 


Let an original high-order linear stable system be represented in the transfer function 
form by : 

to + M 1 + . - + iw’" (4.1.1) 

Go + CLlS 1 + . . . + a n S n 

where a n = 1, a 0 , . . . , a n _ i, 6 0 ? - - * , b m are polynomial coefficients and are probably all 
unknown. 


Let the corresponding reduced model be: 

n ( v _ Co + c^ 1 + . . . + c mo s m ° 

G f ( <s ) — . , 

do + djs 1 + . . . + d no s n ° 

where m 0 < n 0 < n. 


(4.1.2) 


The problem considered here is to set up a criterion for measuring the difference 
between G(s) and G r (s), and to select these coefficients Co, ci, . . . , c mo , do, . . . ,d no 
to minimize the error norm. Here, we utilize the frequency fitting approximation. 
Suppose { (x-’i , u. ! 2 , . . • ,wj} are selected as the fundamental frequency points ranging in 
the system bandwddth. The choice of cq, Ci, . . . , c mo , d 0 , . . . , d nQ should ensure that the 
model frequency response G r (juic){k = 1 , . . . , J) matches the original model frequency 
response G(j&k) as closely as possible. 

One simple approach for fitting the frequency responses at points (uq,u> 2 , . . . ,u >j) 
is to set up the following equation errors: 

C(jwi) - D(jui)G t = e, 

i = 1,2,3,... ,J (4.1.3) 


69 


or equivalently, rewriting the (4.1.3) into the equations of real and imaginary parts. 
Co — C2 Ujf + C 4 ujf + . . . 

—do\G,\ cos {<f>i) + d\\Gi\u>i sin (</>,) + ^ 2 \Gi |u,f cos(<?!>,) + . . . = ef 

ClW, - C 3 U>f + C$U)f - ... 

-<io|Gi| sin(&) - di|Gi|u;, cos(<£.) + d 2 \G,\u:f sin {(pi) -... = £■ (4.1.4) 

i = 1 , 2 , 3 , . . . , J 

where, G{ is a frequency response estimation at u>, obtained through the algorithms 
developed in the last chapter. 

Gi = \Gi\e j<t> ' 

and 

ef = &(£«), 4 = $(£<) 

e, is the frequency fitting error. 


Assume d no = 1, and 

<xf = ^{(M) n o|G,|(cos(^) + j sin(^)} 

o{ = |G,|(cos(^) + j sin(</>,)} 

Collecting all unknowns into the parameter vector: 

Co 

Cl 


£ = 


C m 0 


do 


(4.1.5) 


(4.1.6) 


d 


no — 1 J 


70 



and defining a dimension 2 J x (mo + no + 1) data matrix as follows: 

1 0 -u>j 0 u> j ... -|Gi| sin(cii) |(?i|u;i cos(<^>i) 

— 0 ... — |Gi| sin ( ) — |Gi cos(</>i ) 


0 UJi 


0 


1 0 —u)j 


Lo uij 


0 u* .. 
-^3 0 . . 


— |Crj| sin(4>j) |Gj|u?3cos(^) 
-|GV|sin(<A/) —\Gj\uj cos(^j) 


(4.1.7) 


Define: 



■erf' 







a = 

... ^ 
b 

S — 

$ 




-4- 


(4.1.8) 


Then we can combine the real and imaginary parts of equation(4.1.4) into a stan- 
dard regression format for setting up the least squares estimation: 


rnine = rnin(E£ — a) T (E£ — cr) 


(4.1.9) 


A consideration of this procedure leads to the following remarks: 


1. Estimating the unknown parameter set of the reduced model only involves 
constructing and solving linear equations. Hence, the algorithm is simple and 
the computational cost is low. 

2. The reduced model may have large distortion in some frequency bands com- 
paring with the original system because the selected error norm (4.1.3) differs 
somewhat from fitting the true frequency response. Actually, (4.1.3) can be 
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regarded as the true error criterion weighted by the denominator polynomial 
D(ju). Thus there is a trade-off between the model accuracy and simplicity 
of the algorithm. 

3. The algorithm cannot guarantee the stability of the reduced system although 
the original system is stable. The possibility of getting a reduced stable 
model will be enchanced by picking up more middle and high frequencies in 
the system bandwidth. 


4.2 Scheme 2 

Scheme 1 presents a biased result because it is equivalent to weighting the match 
criterion by the reduced denominator D(ju>). Here, the following iterative algorithm 
can be employed to asymptotically remove the bias if it converges. The error criterion 
at the k th iteration is defined as: 

E t = £ I n W ‘ - . '[ftlM) - G,D„(ju ,) ] I 2 (4.2.1) 

Dk-i{j Wi) 

where D k -i{jui), i = 1,2, . . . , J, are the estimates of D(ju) at the (k- l) th iteration. 
Another similar scheme is to utilize the following error criterion: 

E k = j^W^\e k , | 2 

1 = 1 

= C„(j a.,) - - 11 - ft (4.2.2) 

E k -\ \]Ui) 

These recursive least squares algorithms update the estimates of the reduced sys- 
tem at each iteration. If C k {ju) -* C{j u>) and D k {j uj) D{ju>), the iterative 
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algorithms will asymptotically produce a reduced order system which matches the 
frequency response data without any bias. 


One advantage of using equ.(4.2.1-2) as an error fitting criterion is that we can 
use a linear least squares algorithm to update the estimation at each iteration be- 
cause each error term is linearly related to the unknown parameters. Thus, nonlinear 
optimization algorithms are not needed. 

From the error criterion (4.2. 1-2) we can see that, at the k th stage of iteration the 
error term e*,, has the general form: 

e*, i = Xk,iCk(jui) + pk,iD k {jui) ~ C k,i (4.2.3) 

where for the error norm (4.2.1), 

A Wi 

*’* D k -i(ju>i) 

Wi * 

Pk,i = ~~n T — \ Cj ' 

Di t-i(jwi) 

a, = 0 (4.2.4) 


For the error norm (4.2.2), we have: 


Pk,x = — 


a, = w<[ 


A fc .i = Wi 

Ck-l(jU),) 

Dk-i{jui) 

Ck-l{ju)i) 


Dk-i(jui) 


W, 

- Gi ] 


Define: 


cj = [ CQk c 2 k c 4k . . .] 
C s [ ^3 k ^5 k • • •] 


(4.2.5) 
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d<F = [ dok d 2 k d 4 k ■ ■ ■} 
dj = [ d\k d^k d$k ■ • ■] 

«5 = [10 -^0^0 ...] 
nji = [ o - lj] o - u? o wf . . .] 


(4.2.6) 


Rewrite A*;,,, /£>*,», (k,i, Ck{j>^i), Dk{j^\) as the real and imaginary parts as: 


A k,i = A 


(R) 

k,i 


+ jX 


V) 

k,2 


(R) , • (/) 

Pk,i = Pk,i+JPkj 


Cw = C + KIJ 


dl) 


Ckij^i) — ^c,i c c T 
Dk{j^i) = ^c,i^ c "h 


(4.2.7) 


Therefore, at the k th stage of iteration and at each frequency we can represent 
the error linearly in terms of unknown parameters: 


*(e*,) = X^Cllcc ~ + Pi>l*< ~ ~ Cl? 


9(e w ) = A'^f^cs + A^fiJiCc + P^^ids + - Ckj 


il*) n T 


(«)o T, 




*(/) 


(4.2.8) 


Express the equations (4.2.8) as a matrix form: 


■*£?«» 


T 

c,i 







{Cc > 

-AIM 



T 

r A 

AIM 



d? 




Id?) 


i = 1,2,3, 

...,j 



/d? 

I d? 


(4.2.9) 
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Collecting all these equations, we can construct a standard linear regression form 
for the least squares solution. Without loss of generality we can assume d no = 1. If 
(k — 0, we will get a trivial solution. In this case equation (4.2.9) can be modified 
slightly to overcome this problem. 

The following observation can be made for this scheme: 

1. As in scheme 1, this procedure can be implemented via the linear least squares 
algorithm, thus it is easy to realize. Also it might asymptotically eliminate 
the matching bias if it converges. 

2. Numerical simulation and analysis have shown that sometimes the matrix 
adaptation (4.2.9) will be divergent. Also when noise is not neglectable, it 
will not guarantee convergence to a desired reduced model. 

3. At each iteration, we have to solve an overdetermined set of equations. The 
computation cost is larger comparing to the first scheme. 

In addition, the above algorithm lacks the flexibility of dealing with magnitude 
and phase information. In the following, we develop a procedure matching a reduced 
model and the frequency response data directly. The expense we will pay is by 
employing a nonlinear optimization technique. 


4.3 Scheme 3 

We assume that the simplified model (4.1.2) is specified in terms of second order 
cascade connections and a scalar constant, which is determined by a parameter set P. 
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(G, , G 2 , • • • , Gj) are treated as the original frequency responses at selected frequencies 
(u?i,u> 2 , We choose the r th norm error criterion by: 

E=[t w i IG^-G,# 

t=i 

where, 

Gr,i = Gr (M) 1 = 1,2, . . . , J (4.3.1) 


Our goal is to find optimal parameters P to minimize E, or equivalently E T . For 
utilizing a nonlinear optimization algorithm, we need a formula for calculating the 
gradient. Suppose q € P is a parameter, then it is easy to prove that the derivative 
of E T with respect to q is given by the following equation: 


where, 


rl F T J BG* 

= J>|G r ,< - G t r 2 W^\-^(G r , t - Gt)] 
dq fr{ Oq 


dGU _ dG* r (ju>) 

dq dq 


(4.3.2) 


(4.3.3) 


G* is the complex conjugate of G r . 


For the second order cascade connection form of a reduced model G r , we have the 


following: 


L 


G r (s) - G 

k = 1 


1 ~h 1^ ~f 

1 + Qk, 3 s + qk,4 s2 


CF(s) 


(4.3.4) 


The advantages of choosing the form (4.3.4) are its relative lower parameter sen- 
sitivity and direct relationships between the parameters and the location of zeros as 
well as poles. Based upon (4.3.4), for each parameter q e P we can represent the 
model G r {jijo) as a product of the two terms: 

Gr(ju) = G r ,.{u)g(q,u) (4.3.5) 
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where only g(q,ui) involves the parameter < 7 , G rt s(u>) does not depend on q. Parameter 
q could be associated with the value of a real pole, a real zero, or complex conjugate 
poles and complex conjugate zeros. 


By substituting equation(4.3.5) into (4.3.3), a simple calculation will lead to the 
following equation: 


where, 


= X>|G r , t - " Gir 3 W&[G* r j(G r ,i - G,) aing -~ ] 
oq oq 


q.iq) = g(q,<*>) 


(4.3.6) 


(4.3.7) 


If we select r to be 2, and take the derivative of E 2 with respect to the scalar pa- 
rameter C and set the equation to zero, then the parameter C can be easily estimated 
in terms of the other remaining parameters. 


In fact, using (4.3.2) we can obtain: 
f)F 2 J J 

— = 2 x; - G,)] = 2 £ Wi[C | F, I 2 -R(F*G,)] = 0 

t=l «=1 


(4.3.8) 


where, 


Fi = F(ju) 1 ^ 


Solving the equation for C, we get: 

c ZL WiX(F?Gi) 

Ei, Wi I F I 2 

In this way, we can reduce by one the total number of unknown parameters. 


(4.3.9) 


Another alternative is to represent the reduced model by taking the partial fraction 
expansion. For each parameter q £ P we can represent the model G r (s) as a sum of 
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two terms: 


G r (ju) = G r ,»(w) + 9i{ u> )92{q, u) (4.3.10) 

where G r , a { w) and ^(w) do not depend on the parameter g, only g 2 {q, w) involves 
the <7. Similar to the previous case, q can associate with a real pole, a real zero, or 
complex conjugate poles and complex conjugate zeros. g 2 (q,u>) could be expressed 
explicitly for all these different cases. 


It is clear to see that the error norm criterion defined in equ.(4.3.1) takes care of 
both magnitude and phase information of the original system. If only the magnitude 
is important, we can further simplify the model structure and keep the quality of 
matching magnitude response as good as the previous one. The error norm criterion 
we can use is: 

E = EH'illGr.il - l&iri’ (4.3.11) 

t = l 


It is easy to show in this case that the formula for calculating the differentiation 
is given as follows: 

dE ’ - |g.ll r ' 1 |G,,^( a '°^ - ( --)] (4.3.12) 


dq 


1 = 1 


here, all the variables have the same meaning as the former equations. 


It is well understood that it is not guaranteed for a nonlinear optimization scheme 
to achieve the global minimum in the general case. Actually, it is most likely for 
the algorithm to get into a local minimum if the initial value is arbitrarily selected. 
Therefore, in order to get a good estimate of the optimal reduced model, it is necessary 
to have a fairly accurate initial estimate . In the next section, we will present a quite 
reasonable algorithm to estimate the initial value for the reduced system G r (jui) which 
will help the optimal algorithm obtain excellent results. 
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The algorithm shown here could be utilized to weight different frequency bands, 
and force the system to satisfy a stability constraint. Usually, the deviation |G ri , — G , \ 
or — | G,'|| differs siginificantly in different frequency ranges. The error at the 

middle frequency band tends to be larger. Hence, utilizing adaptive weighting factors 
W t is preferred instead of using constant values. Corresponding to the criterion 
(4.3.1), the usual selection of adaptive weighting factor is: 


w,= 


|G r ,j — G,| 

£i, I G„ - G, | 


(4.3.13) 


For the error norm (4.3.11), we choose: 

HG.31 - fen 
‘ Eh lia.ii-iG.il 


(4.3.14) 


During the search, W x will be modified adaptively according to (4.3.13-14) with the 
updating of all the parameters q £ P and G> ?t in each iteration. 


It should be mentioned that the selection of the fitting frequencies is very im- 
portant. Any inapproprite selection of dominate fitting frequency points will lead 
to a bad accuracy in some frequency bands which might be crucial to the reduced 
model. One thing obvious is that the selected frequency points should cover the 
system bandwidth and the total number of fitting frequencies Lj should satisfy : 

Lf > n 0 + m 0 + 1 

Assuming that the estimated frequency response G : is exact, a larger Lj implies 
the reduced model will be more accurate. The expense we pay is to increase the 
computation cost. 
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4.4 An Algorithm for Initial Value Estimation 


In this section, a least squares algorithm is suggested to estimate the initial value 
(c 0 , ci , . . . , c mo , d 0 , . . . , d no — i) of the reduced system. 

Construct the following equation error based upon the complex sinusoidal modu- 
lating signals: 

m+ n 

E3 b k - m [D(-jku 0 )Y k - C(-jku) 0 )U k ] = £ m (4-4.1) 

fc=m 

where, 

C( — jkuJo) — Cq + Ci ( — j kid o) + C2( jkuj o) + . . . + c mo ( jkldo] 

D(—jku 0 ) = d 0 + d\(—jku!o) + d 2 ( — jku> 0 ) 2 + • • • + d no (—jku { 0 ) n ° (4-4.2) 


Rearranging equation (4.4.1 ) into real and imaginary parts, we get the forms: 


m + n m+n m+n 

do X^ bk- m^k + ^1 Xw h-m(kuJ 0 )Y k S — C ?2 ^-771(^0) Yk + • * • 

k=m k—m k—m 


m+n m+n d 

—Co E^ b k - m u k — Ci E^ b k - m (kido)U k + c 2 EZ b k ^ m (ku>o) 2 U k + . . . = e m 

/c=m k=m k=m 

m+n m+n m+n 

do £ 6*- m n s - di Y2 h- m {kdJ 0 )Y k c b k . m {ku 0 ) 2 Y k s + . . . 

k—m k=m k—m 

m+n m+n m+n 

— Co EZ b k -mU k + Cl Y2 b k -m(kuJo)U k + C 2 EE b k - m (kijJ 0 ) 2 U k + . . . — £ m 
k=m k—m k=m 


m = 0, 1, 2, 3, . . . , M 


(4.4.3) 
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oe o <c? 


Define: 


0v = 


r d 0 

d i 

dno— 1 
Co 
Cl 

c mo J 


For each time interval + T], we assume that: 


(0 = 


■ zufemo ■ 


' ELoM *-*)W) - 

HUWW) 


-ELo Mfc*)»?(*) 

ej±! H-in'(o 

r?(«) = 

ES 6*., (*■*>)!?(•) 

ees 



: 

2(M+1) 

: 


J 2(M+1) 


(0 = 


'* 1+1 
^k = 1 

>n+l 


WHO • 


• -ZUoh(kw 0 )UHi) ■ 

MW) 


+ Et=o hih^WUi) 

*-.lW) 

rj(0 = 


*-.tW) 

na/ i i \ 

+ EIt. , M.(*u’B)tW) 


J 2(M+1) 


2(M+1) 


Define: 


r(») = [r S(0, r?(z), . . . , r^, r“(?:), . . . , r“ (i)] 


Assume: 

m-j-n 

*£ = *< z + jy;(,)) 

k=m 


(4.4.4) 


(4.4.5) 
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(4.4.6) 


m+n 

/i = i l ..H^r(n , (i)+ in*(0) 

k=m 

m = 0, 1,2,3, ...,M 


and denote: 


' Mo (*) ‘ 

Mo(0 


C(0 = 


Mm(0 

Lmm(0 


2(M+1) 


(4.4.7) 


Collecting all the equations for m = 0, 1 , 2, . . . , M, we form the following matrix 
equation for each time interval [f,-, C 4- 7 j : 

e(t) = r(i)0 P - C(0 

f = 1,2,3,..., W (4.4.8) 


Make the following definition: 



■ r(l) ' 


■ C(1) • 

r = 

r(2) 

C = 

<(2) 


S(N). 


-C(A-). 


Collecting all matrix equations (4.4.8) for each time interval [t,, t;+T], we construct 
a standard linear regression equation for the least squares estimation: 

mine = min(r0 p — C) T (r0 p — () (4.4.10) 

e p e„ 


Equation (4.4.10) can be used to estimate an initial value for a reduced system. 
Again there is no promise for these initial values to satisfy the stability constraint. 
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The initial reduced system matches the harmonic Fourier series coefficients of the 
input and output data. The estimated initial values are fairly accuracy. These values 
will be updated to satisfy the optimal criterion. 


4.5 Some Numerical Examples 


In this section we present some numerical examples to demonstrate the performance 
of the algorithm by comparing with the existing approaches. Here we apply the 
scheme 3 for model reduction. For the MFFM algorithm, there is no requirement of 
knowing the actual high order system a priori. But we do need an input-output pair 
generated from the original model. These data can be obtained from the real process. 
If the original model is known, these data can also be produced by simulating the 
system. 

Example 1. Co and Ydsti[CoTB90] compared a Frequency Fitting Pade algorith- 

m and the P-L Modulating Function method in reducing some high order models into 

lower order systems. For easy comparison, we took the first example from [Xihe87]: 

(5 -f 2) 2 (s + 5) 2 (s + 1000) 
ll5J (s + l) 2 (s + 10) 2 (s + 100) 2 

The PL method which matches the Fourier series coefficients of input-output data 

is quite similar to our scheme 1. An input signal for exciting the system used by 

[CoTB90] consists of a positive unit step function followed by a negative step function. 

Data was collected in five seconds. The reduced third order model is: 

PL _ 3. Sills 2 + 21.78535 + 38.7277 

1 ^ “ s 3 + 49.9836s 2 + 548.5900s + 381.0047 

They found that the PL method works better than the FFP algorithm. 
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For the purpose of comparison, we also utilized the classical balanced realization 
method to reduce the system. Using functions provided by MATLAB, we obtain a 
reduced model: 

Rln 5.4156s 2 + 30.9466s + 58.5745 

Hl ^ “ s 3 + 65.4461s 2 + 793.1840s + 585.7454 

For simulating our MFFM algorithm, we selected the same random binary input 
signal as used in chapter 3 to excite the system. Data is collected in 8 seconds. The 
whole data set is divided into the six time intervals [<i, ti + 3] for t x = 0, 1 , . . . , 5. The 
sampling time interval At is 0.005. 1024-FFT is utilized in calculating the Fourier 
series coefficients on each time interval. We first estimated the frequency response 
at frequency knots { kuJo , k — 0, 1,2, . . . ,30}, uio = 2.094. 

The reduced third order system by matching the first twelve frequency knots is: 

T r MFFM 4.5233 s 2 + 30.6739s + 45.2871 

H 1 _ s 3 + 51.2456s 2 + 710.2395s + 453.6771 

Fig.4.1(a)(b) present the Bode plots of above system along with the simulation results 

obtained by applying the balanced realization method and PL procedure. Fig. 4. 1(c) 

shows the corresponding Nichols plot. 

The low order model achieved by matching the last fifteen frequency knots is: 

- , 1.4952s 2 + 863.1844s + 5334.4197 

1 1 ; ” s 3 + 187.4456s 2 + 11304.2301s + 141250.0129 

The Bode plots of the above system along with other results are shown in Fig.4.2(a)(b). 

Fig. 4. 2(c) gives the corresponding Nichols plot. 
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Frequency(rad/sec.) 


Fig. 4. 1(a) Magnitude response vs. frequency. The MFFM uses the low frequency 


matching. 



Frequency(rad/sec.) 

Fig. 4. 1(b) Phase response vs. frquency. The MFFM uses the low frequency 


matching. 
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Magnitude Response(dB) 



Phase(degree) 


Fig. 4. 1(c) Nichols plots. The MFFM uses the low frequency matching. 



Frequency(rad/ sec.) 


Fig. 4. 2(a) Magnitude response vs. frequency. The MFFM uses the high frequen- 
cy matching. 
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Frequency(rad/sec.) 


Fig. 4. 2(b) Phase response vs. frequency. The MFFM uses the high frequency 
matching. 



Fig. 4. 2(c) Nichols plots. The MFFM uses the low frequency matching. 
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Simulation results indicate that the PL method and the balanced realization tech- 
nique weight the low frequency modes more heavily. The MFFM algorithm gives 
results similar to the PL and balanced realization by picking more low frequency 
knots. If we sacrifice the DC response and the frequency band lower than 10(Hz), 
then the MFFM alogrithm gives a reduced model which approximates the original 
system much better in the frequency band larger than 10(Hz). The tradeoff depends 
upon real applications. The MFFM has this flexibility. 

Example 2. The second example we consided here is a six order model: 

4.5000s 5 + 16.8750s 4 + 1.2474 x lOV + 8.0454 x 10 6 s + 3.6556 x 10 7 

= s 6 + 66.85s 5 + 2778.9s 4 + 7.1963 x 10 4 s 3 + 1.0168 x 10 6 s 2 + 9.8011 x 10 6 s + 4.7306 x 10 7 

The same random binary signal is used to excite the system as in Example 1. 
Data was collected in 8 seconds. The w r hole data set is divided into two disjoint 
time intervals for using the PL reduction. The parameter L was selected as 12. The 
reduced system is: 

PL 0.1372s 2 + 124.0337s + 345.1078 

~ s 3 + 10.6977s 2 + 145.0830s + 453.2253 

Llsing the balanced realization, we obtained: 

o, -5.5419s 2 + 170.4683s + 872.5521 

Tjtfln/ \ ___ — — 

2 [ ’ s 3 + 15.4064s 2 + 207.4581s + 1129.1082 

For using the MFFM method, we divided the whole data set into the six time 
intervals [t„ t, + 3], f, = 0, 1, . . . , 5. We first estimated the frequency response H 2 (j u?) 
at knots { ku 0 , k = 0, 1, 2, . . . , 30}, u 0 = 2.094. We picked five low frequencies, five 
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middle frequencies near the valley and five high frequency knots. We put a slightly 
heavier weight on the middle frequencies. The reduced third order system is: 

jjMFfm M 3.8121s 2 - 59.6845s + 6913.6012 
2 _ s 3 + 42.7534s 2 + 389.3838s + 8388.2031 

Fig.4.3(a)(b) present the Bode plots of obtained results. Fig. 4. 3(c) is the correspond- 
ing Nichols plot. 
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Frequency(rad/ sec.) 

Fig. 4. 3(a) Magnitude response vs. frequency. 



Frequency(rad/sec.) 

Fig. 4. 3(b) Phase response vs. frequency. 
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The MFFM achieves considerably better results in this example. The PL algorithm 
produces a large bias in the high frequency band. The balanced realization cannot 
track the valley. Only the MFFM gives almost an exact match for low frequency and 
high frequency modes, as well as tracks the rapid change in the frequency band near 
the valley. 



Fig. 4. 3(c) Nichols plots. 
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Chapter 5 


High Resolution Frequency Estimation in 
the Presence of Noise Using Complex 
Sinusoidal Modulating Signals 


5.1 Introduction 


In this chapter, we focus on the problem of estimating the angular frequencies u>i,u? 2 , 
. . . ,l:k for given superimposed harmonic signals over a set of finite time intervals 

{ [tiiti + T], i — 1, 2, . . . , Q) ■ 

y(t) = Vd(t) + e(t) (5.1.1) 


where, 


K 

yd{t ) = Y. A ' sin ( u; «^ + 

t =i 


(5.1.2) 


In equations(5. 1.1-2), e(<) is a stationary white Gaussian noise, K is the number 
of superimposed sinusoids and is not necessarily known. Ai and <f>{ are unknown am- 
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plitudes and phases of the signal yd(t). T, the length of the time intervals, is actually 
a frequency resolution related parameter. The time intervals need not necessarily be 
disjoint. 

The problem has many practical applications and has received considerable atten- 
tion[KayS84, 0hSG91, Sher91]. Although traditional FFT and periodogram methods 
work well in a high signal-to-noise ratio case, they are nonparametric algorithms and 
the resolution is not high enough for small data sets, nor for low SNR’s. 

This has spurred interest in developing high resolution parametric algorithms. 
Some approaches, like the maximum likelihood estimation, nonlinear optimization al- 
gorithms, etc., work well in the low SNR situation[Fraz88, Scha91, Tuft82], However, 
they require either long data lengths or sophisticated computing and time costly algo- 
rithms. MUSIC is a typical algorithm for frequency estimation and direction finding. 
But it is quite computationally demanding[Stoi91, LeeH91]. The Yule-Walker equa- 
tion approach has been widely utilized for frequency estimation and spectrum analysis 
because of its effectiveness and the availibality of fast lattice algorithms. Recently, 
some researchers have proposed a HOYW algorithm for improving the frequency res- 
olution [Chan82, Mose88, Stoi91] . But the problem is that sometimes it is difficult 
to separate “ spurious zeros ” from signal frequencies. 

In contrast to the previous methods, here we focus on the problem of high reso- 
lution frequency estimation employing a simple linear least squares algorithm. We 
use a continuous time model for the analysis. The method we utilize is to construct 
an autoregressive differential equation model to fit the received signal. Then com- 
plex sinusoidal modulating signals are utilized to convert the differential equation 
into simple algebraic equations. The parametric least squares algorithm can be de- 
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signed to detect the frequencies. One advantage of our method is that there is no 
requirement of a long data length record, i.e., the algorithm works well for short time 
observation data. Numerical simulation results will demonstrate the performance of 
the algorithm. 


5.2 Overview of the High Order Yule- Walker 
Estimation 


The HOYW estimation utilizes the discrete time samples of the sinusoidal signal 
(5.1.1) for t = 0, 1, . . . ,N — 1. Assume ^ and the phases <j) ,, i = 1, 2, . . . , A 
are mutually independent random variables obeying a uniform distribution on [0, 27r], 
with frequency uik £ [0, 2 7 r ] . It is assumed that the signal j/d(t) and the white Gaussian 
noise e(f) are uncorrelated for all t,L i.e., E{yd{t)e(t)} = 0. 

Define a whitening filter as: 

A{z~ l ) = IK 1 - z'^Kl - 

fc=i 

= El ( 1 — 2z -1 cos(uijt) + z~ 2 ) (5.2.1) 

k=\ 

Then it is easy to verify that the received sinusoidal signal y(t) obeys the following 
discrete ARMA model: 

A{z-')y(t) = A(z~ l )e{t) (5.2.2) 

The filter has zeros on the unit circle e ±J ^, k = 1,2, . . . , K, from which we can get 
the frequency estimates of the signal once the zeros of A{z~ 1 ) have been estimated. 

It is easy to verify from equation (5.2.2) that the signal y(t) also satisfies the 
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following high order ARMA model: 


C( 2 -')y(t) = C(z-')e(t) 


where, 


and where B(z *) 


C(z 1 ) = H C * 2 ‘ = B( z X )A(z *) 

*=0 

is a polynomial of degree L — 2K\ 

L-2K 

B(z-')= E 

t=0 


(5.2.3a) 


(5.2.36) 


(5.2.3c) 


It has been shown from theoretical analyses and numerical experiments that using 
the high order ARMA equation (5.3.3a) for frequency estimation can decrease the 
variance of the estimator and improve the frequency resolution. But the problem is 
that sometimes it is not easy to separate the “ spurious zeros ” of B(z~ l ) from the 
signal frequency zeros of A(z -1 ). 


Define {r(fc)} to be the autocorrelation sequence of the signal, i.e., 


r(k) = £{»(/)»« + *)} 


(5.2.4) 


Without loss of generality, cq can be constrained to be 1. It is well known that the 
following high order Yule- Walker equation holds for solving the coefficients of the 
whitening filter: 

/ r(L + 1) r(L) ... r(l) \ / Co \ 

r(L + 2) r(L+ 1) ... r(2) 


c 1 


= 0, L, Li > 2K (5.2.5) 


\r(I + I,) r{L + h- 1) ... r(L\) / \c L f 


The number of equations L\, and the order L of the polynomial C{z) can be 
selected flexibly for improving the resolution and accuracy of the frequency estimates. 
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But the larger L is to be chosen, the more difficult it is to separate the roots of B(z) 
and A(z). 


Rewrite (5.2.5) into the following standard linear regression equation(co = 1): 

(ci c 2 ... Cn) R(L,Li) = - (r(L + 1) r(L + 2) ...r(L + Li)) (5.2.6) 


w T here, 


(r{L) 


R{L,L ,) = 


\r(l) 


r(L + Li — 1) \ 

r(-^i) / 


Numerical experiments have shown that the covariance matrix R{L.L\) is very 
likely to be ill conditioned, especially for large (L,Ii). For achieving a good robust 
estimator, some regularization procedures have to be taken. It has been proven that 
a much better frequency estimate can be obtained by first approximating R(L , L\) in 
the subspace of rank 2I\ in the Frobenius norm sense, and then taking the Moore- 
Penrose pseudoinverse. See [Stoi91, Chan82] for more details. 


The HOYW equation has some advantages. The efficient lattice Levison algorithm 
can be utilized to solve equation(5.2.5), and computation time can be saved. The 
selection of (L,Li) is flexible. In the next section we present a method of using 
complex sinusoidal modulating signals to estimate the frequencies. Comparing with 
the YW equation method, our approach gives high frequency resolution estimates, 
especially for short data lengths. The proposed algorithm is very robust and shows 
outstanding performance in high SNR circumstances. Conversely, it also works well 
in the low SNR case. 
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5.3 The Parametric Least Squares Estimation 


For the received signal (5. 1.1-2), we assume e(t) is a real stationary Gaussian white 
noise with zero mean, and variance a 2 . Our estimation algorithm is based upon the 
fact that the signal yd{t) = sin(o>,t + <f>) satisfies a differential continuous time 
autoregressive model of order 2 K. The coefficients of this model only depend on 
the angular frequencies and not on the amplitudes nor phases. More specifically, K 
imaginary zeros of the autoregressive model are exactly the K angular frequencies to 
be estimated. 


Define: 


n(p) = n£ =1 (p 2 + u,£) = p 2K + Qjp 2 ^'- 1 ) + . . . + a A _ lP 2 + a K 


where, p = d/dt(differential operator), 


"i = 

J = 1 

Q 2 = u i u k 

J*k 


(5.3.1) 


K 

ctK = n S 2 (5.3.2) 

j=i 

Qi, a 2 , . . . , Q k are elementary symmetric functions of loi,u) 2 , . . . ,ujk. 


Then it follows that on each time interval [t,, t t +T], the received signal y(t) satisfies 
the differential equation model: 

n (p)y{t) = y 2h + Ojj/ 2 ^ -1 ) + . . . + a h -_ iy 2 + a K = e(t) (5.3.3) 

where, 

t{t) = n(p)e(t) 
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and e(t) is the noise term. From (5.3.1) we can see that U{p) has imaginary zeros at 
tn — 1, 2, . . . , Fl . 

Define a set of modulating functions: 

<f> m (t) = e jWot (l - e JUJot ) n (n > 2 K) (5.3.4) 

m = mo, m .\ , m 2 , m ; v/ 

where u> 0 = 27 t/T. The appropriate selection of these m, will make the modulating 
signals more flexible and work in different situations. Taking the binomial expansion 
of (5.3.4) and changing the index of summation, <p m {1) can be written as the form: 

k - 0 

6 t = (-l)‘Cj (5.3,5) 

where C k is the binomial coefficient. 

Multiplying both sides of equation(5.3.3) by <j> m (t) and using the modulating prop- 
erty, we can construct the following equation format for the observed signal on each 
time interval -f T], t = 1, 2, . . . , Q: 

m+n 

E b k - m n{-jku 0 )Y k {i) = e m (r) (5.3.6) 

k—m 

m = m 0 , mj, . . .,m M . 

where 

n(*)= [ T y(ti + t)e jk “ o1 dt 

Jo 

m + n 

c m (0 = E 

A:=m 

£*(*)= [ T e{t t + t)<j> m (t )di (5.3.7) 

JO 
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The best selection of m 0 , mi, , wm is to make [rriju> o, (m ; + n)u,’o] cover, or at least 
be close to, the signal frequencies. In this way eqn.(5.3.6) will most reveal and utilize 
the signal information. 

From equation (5.3.1), we have: 

n(— j kuJo) = qk — (fcwo) 2 OTi'-i + • • • + (jku> o) 2 ^ J) ai + {jku> o) 2h (5.3.8) 

and hence rearranging equation (5.3.6) into a linear form relating the unknown pa- 
rameters a,, we obtain: 

m + n m+n 

CtK Y h-mYk(i ) - Qft'-l Y h-m(kuo) 2 Yk{i) + ■■ ■ 

k=m k=m 

m+n m + n 

+ (-l) A - 1 a 1 Y b k . m (ku 0 )^Y k (r) + Y h- m (-l) K (ku> 0 ) 2K Y k (i) = «m(») 

k=m k=m 

m = mo, mi,m 2 , . • ■ , mm 

* = l,2,...,g (5.3.9) 


Define: 


/ <*K \ 
&K - 1 


Q = 


m+n m+n 

Emit) = [ Y h-mY k (l), - Y h-m{kuofY k {i) . . .] 

k=m k=m 

m+n / Y£(i) \ 

/> m (0 = (-i) A 'E w^o) 2A ' 

k=m \Y k (i)J 


(5.3.10) 


Collecting all equations for m = m 0 , m\, . . . and constructing a standard 
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regression matrix-vector format for least squares estimation: 


f =o(0 \ 
Si(0 

OL + 

( Po{i) \ 
p i(0 

\ Hjif(i) / 


\p\f(i)f 


* = 1,2,. ..,Q 


(5.3.11) 


where 


( «o(*) \ 
*i(0 


Using standard least squares estimation, we can get estimates of the coefficients 


di.Q 2 , • • • for the autoregression model. Then we can easily obtain the frequency 


estimates using the following procedure: Construct the polynomial equation. 


X 1 ^ — Q\X k 1 + 2 — ... + (— 1 )* ok — 0 


(5.3.12) 


and solve this equation to get K roots, say xi, £ 2 , • • • 1 3 'a- Then we obtain the fre- 
quency estimates by taking the square root of these roots: Ui = -y/x, (* — 1, 2, . . . , A ). 


The polynomial (5.3.12) is called the Prony polynomial in numerical analysis. 
The accuracy of the frequency estimates will depend upon the estimates of the model 
coefficient cb and the root estimates of equation (5.3.12). The former is more essential. 
If any roots of (5.3.12) are negative or complex, they are regarded as outliers. In this 
case, the whole procedure fails to give the good results. In the next section we prove 
that the variance of the estimation error goes to zero as the signal-to-noise ratio goes 
to infinity. 
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5.4 Performance Analysis for Large SNR 


Circumstances 

In the last section, we have proposed a parametric least squares algorithm for esti- 
mating frequencies of a received signal. Here, we will give the performance analysis 
for the parametric least squares method in large SNR cases. For large SNR assump- 
tion, we mean Ak >> o for each k. So, the noise term e(t) in the received signal 
y(t) = y<i{t) + e(f) is a very small component. Also, we make an assumption here 
that there is no numerical calculation error involved. Estimation bias is only caused 
by noise. In our analysis, we will approximate the original autoregression model and 
discard some high order small terms relating to noise e(t) whenever appropriate. 

If there is no noise, i.e., e(t) = 0, the received signal y(t) will obey the au- 
toregressive model n(p) exactly. When a noise e(t) exists, but the SNR is very 
large, the linear regression equation will give a coefficient estimate o, and a — a 
is a small term due to the least squares estimation property. Solving the roots x 
of the polynomial equation (5.3.12) and taking the squares root of x will give the 
frequency estimates u>k,k = 1,2,...,A\ The estimates will involve a small error 
Sk = u?l — xic, k = 1,2, . . . , K. Thus in the general case we can say that the estima- 
tion errors 6k are random variables depending upon the noise and are small terms 
comparing to the signal. 

In the following analysis we define: 

K K 

n(p) = II (P 2 +U 2 m)= n(p 2 +^m+^m) (5.4.1) 

m= 1 m=l 

Based upon the above argument, we can claim that the following equation approx- 
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imately holds: 


(5.4.2) 


n(p)[y(0] - 0 


Taking the Taylor expension of (5.4.1) and neglecting the high order small terms, 
we have: 


K 


(5.4.3) 


where, 


n(p) = n(p) + 52<5 fc n fc (p) 

k = 1 


n*(p) = flip 2 +“J) 

i=i 


Discarding the high order small terms relating to e and 6 k , we obtain the following 
equation: 

n (p)[y(<)] = E ^njfc(p)[yJ + n(p)[e] = o (5.4.4) 


Using the parametric least squares method, we multiply both sides of (5.4.4) by 
the modulating function and employ the modulating property: 

K m + n 


A m + n r 

EtE bk-mni{-jkwo)Pk]6k + J n[e(t)]<t>m{t) = 0 


(5.4.5) 


/= 1 k—m 


where 


Define: 


ft k — f E Aj sin(w,t + <fit)e 3k “ ot dt = E 
■'o ,=i »=i 


n,(-jfcw 0 ) = n* =1 (a ; 2 - fc 2 u; 2 ) = f I/ W 

T^l 


(5.4.6) 

(5.4.7) 


E k = [ e(t)e jku,ot dt 
Jo 

We can see that the E k are complex Gaussian random variables with E(E k ) = 0 
Var(£’fc) = cr 2 T. 
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Assume: 


m + n 

~ ^ ^ bk~m n ( J kuJQ ) Efc 

k—m 

Utilizing these notations, equation (5.4.5) can be expressed by the format: 

K m+n 
/= 1 k=m 

m = mo,mi,m,2, . . . ,ttim (5.4.8) 

We can represent the equation(5.4.8) as a standard linear vector-matrix form. If we 
choose m, = i, then the regression form is very simple. 


Define: 


B = 


and assume: 


r = f 


def 
* = 


Mi ••• 

bnPn 


\ 

Mi 

^n— 1 

bn $n + 1 



bopM 

... J 

Vo.i 

Vo, 2 

7/O.A- \ 


Via 

Vl,2 

J/l.A 


V VM+nA VM+n,2 

- V-M+n,K ) 



Eq \ 




E i 



€ = 





V -EW+n / 



bi*i ... 



\ 

bo'i'i ... f 

*n— 1 

bn ^n+l 



Mm 

• • - ^n^M+rj / 


(5.4.9) 
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where 


'I'm = n(-jmwo) 


m 


0,1,2,. ..,M 


Thus, equation (5.4.8) can be written as a standard matrix-vector regression form: 

Br6 = (5.4.10) 


/ 6 1 \ 

*2 

V 6 k ) 

It is easy to prove that B has full rank. So if r has full rank, we can solve for 6 
explicitly: 

6 = -[(5r) T (Br)]- 1 (5q^)f (5.4.11) 

$ is a frequency estimation error and Equation (5.4.11) reveals an approximate 
relationship between the estimation error and the noise. It is easy to see that E(6) - 

0. Thus our estimation algorithm gives unbiased estimates. 

Also from above analysis we can see that r is a constant matrix. The elements of 
B depend upon the amplitudes A t . The variance of e is 0(c r 2 ). Thus, for large SNR, 

1. e., Afc/cr — » ► oc for k = 1, 2, . . . , A', the variance of 6 will approach the zero matrix, 

i.e., Var(6) 0. 

The analysis up to now demonstrates that the frequency estimation errors are 
random variables depending upon the noise e. The SNR will determine the variance 


Here, 


6 = 
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of 8 . When the SNR is very large and approaches infinity, the variance of 6 is small 
and goes to zero. 


5.5 Some Numerical Examples 


Here we present some numerical examples which compare the modulating signal ap- 
proach with the HOYW equation method, and evaluate the utility of the proposed 
algorithm for improving frequency estimation in various SNR circumstances. 

Example . In this example, the analog angular frequencies, amplitudes and pha 
ses are given by: uj aj i = 50,u> ai 2 = 55,u; ai3 = 100; A\ — A 2 — A 3 = 6; <j>\ — 0. 123. 
c p 2 — 0.541, <ps = 0. The sampling frequency f s is 100(Hz), measuring time T = 0.64 
seconds. So we have a total of N = 64 points. The normalized frequencies are 
fi = 0.079 Hz, f 2 — 0.0875// / 3 = 0.1592 Hz respectively. The represent the 
normalized angular frequencies, i.e., u^i = 0.5, = 0.55, u; 3 = 1.0. For testing 

and comparing the performance of the algorithms, we select the noise variance a 2 
for producing the different signal-to-noise ratios. SNR t is defined for each frequency 
component /, as SNR t = 10 log 10 (,4 t 2 /2<j 2 ). 

K is chosen to be 3. u? 0 = 27 r/T = 9.82. We use one short time estimation and 
make the very common selection m t = i . The integration(3.7) is approximated by 
the standard parabolic rule. Each m, corresponds to a complex equation. Real and 
imaginary equations basically reveal similar frequency information of the data. So 
we can think them as one equation. The number of unknowns is 3. The minimum M 
needed is 3. The selection of M is important. Based upon the numerical experiments 
and our experience in using the modulating signal approach, we have found that the 
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optimal choice of M is around 6, i.e., making the number of equations double the 
total unknown parameters. This can be explained by the fact that as M increases 
from the minimum value needed for solving the equations, more data information is 
utilized until up to a certain frequency value beyond which frequency more noise will 
be involved in the data causing a deterioration of the LS estimator, see[Pear85]. The 
algorithm gives unbiased estimates in the case of SNR = +oo. 

The normalized error criterion for the frequency estimates is defined as: 
VOEj(dB) = 101o glo (f - *ij) 2 ) 

i=l 

where N is the total number of sampling points, j indicates the j th frequency, u> 3 
is the true frequency and ujij is the estimate for the i ih Monte-Carlo run. Iq is the 
total number of runs. Based upon the theoretical analysis, the asymptotic normal- 
ized Cramer-Rao lower bound(CRLB) is governed by the equation (24<r )/(A?JV 2 ), 
see [Stoi91], The CRLB can be utilized to test the performance of the proposed 
algorithms. 


Tables 5.1 and 5.2 show the normalized variance of d> 3 estimates for several 
selected values of M and SNR based upon the total of 50 Monte-Carlo runs. From 
the tables we can see that the optimal selection of M is 6. As M becomes larger than 
9, the VOEi increases dramatically. For large SNR cases, the sensitivity in choosing 
the different M is low. From the tables, we can see that there is less than a 10(dB) 
difference in picking M — 6 and M = 9 for the SNR— 30(dB) or SNR=40(dB) cases. 
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TABLE 5.1 


NORMALIZED VARIANCE OF C) X (IN DECIBELS) 


M 

STVRj = 0(dB) 

10(dB) 

20(dB) 

30(dB) 

40(dB) 

3 

-5.7335 

-7.9588 

-21.9382 

-28.4797 

-49.9441 

6 

-11.4818 

-23.3867 

-27.5682 

-35.8743 

-54.6284 

9 

-11.3892 

-15.9320 

-19.3961 

-30.4661 

-48.8173 

12 

-9.8970 

-13.9034 

-16.6418 

-27.7312 

-46.661 

14 

-6.3752 

-12.7231 

-16.0842 

-27.6992 

-45.551 


TABLE 5.2 

NORMALIZED VARIANCE OF (IN DECIBELS) 


M 

SN^ = 0 (dB) 

10(dB) 

20(dB) 

30(dB) 

40(dB) 

3 

7.6042 

4.9103 

2.5681 

-6.3952 

-10.7561 

6 

-11.3122 

-26.1734 

-32.8596 

-38.6037 

-57.4839 

9 

4.0824 

-12.3958 

-20.6769 

-30.2958 

-48.5842 

12 

10.5268 

7.4328 

-14.5609 

-20.5128 

-39.4813 

14 

13.6248 

12.9028 

9.5663 

-3.9128 

-21.4812 


Regarding the HOYW algorithm, the fact that we only have 64 points and 3 
unknown frequencies will lead to a poor autocorrelation estimate and increase the 
difficulty in attaining frequency root separation if L and L\ are too large. We selected 
L = 18, L\ = 20 in using the HOYW equations. This is a fairly good selection in 
this case. Fig. 5 . 1 (a-c) and Fig.5.2(a-c) present the normalized variance and bias of 
u>i, Cj 2 , W 3 , using the proposed modulating signal approach and HOYW equations 
based upon the total of 50 Monte-Carlo runs. The theoretical Cramer-Rao lower 
bound is also shown for easy comparision. From the simulations we can see that 
our algorithm is much more robust and gives reasonable results even for negative 
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SNR(dB). Comparing to the HOYW equation method, the variance of the frequency 
estimates of our modulating signal approach is much close to the CRLB in most cases. 
For the case of SNR=-10(dB), approximately one out of eight runs failed because of 
complex or negative roots. These were regarded as outliers and were removed from 
the least squares estimation. For SNR=— 5(dB), roughly one out of fifteen runs failed 
to give a good result. For SNR=0(dB), only one run broke up. 
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Variance of Estimation(dB) Variance of Estimation(dB) Variance of Estimation (dB) 


10 


i i 



ig.5.1(a) The variance of estimation of u?i in different SNR cases. 



Fig. 5. 1(b) The variance of estimation of u >2 in different SNR cases. 



Fig. 5. 1(c) The variance of estimation of l>s in different SNR cases. 
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20 



SNR(dB) 

Fig. 5. 2(a) The estimation bias(dB) verses different SNR for u>\. Solid line 
for the modulating function method, dot line for the HOYW. 



Fig. 5. 2(b) The estimation bias(dB) verses different SNR for u> 2 . Solid line 
for the modulating function method, dot line for the HOYW. 



Fig. 5.2(c) The estimation bias(dB) verses different SNR for w 3 . Solid line 
for the modulating function method, dot line for the HOYW 7 . 
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The HOYW equation method cannot identify the normalized frequencies 0.5 and 
0.55 for SNR smaller than 5(dB). Fig. 5. 3 demonstrates the power spectrum of the 
identified whitening filter(mean value) for SNR = 0(dB). From the graph we can see 
that there is only one peak around the middle of the 0.5 and 0.55 frequencies. The 
HOYW algorithm mixes the two frequencies and only gives one frequency estima- 
tion. Also, there are several spurious frequencies having larger power than the true 
frequencies, which makes the root separation of A(z) and B(z) much more difficult 
and even impossible. 



Fig. 5. 3. Power spectrum of the estimated filter 1 /C(z *) using the HOYW 
equations for SNR=0(dB). 
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If the FFT or periodogram is utilized to identify the frequencies, the minimum 
resolving frequency is guided by the formula A / = fs / ^’ • which is 1.56(Hz) or 
9.8(rad./sec.) in our case. So the frequencies 50 and 55 will be recognized as one 
frequency and cannot be distinguished. 

From our simulations, we have also found that the selection of w 0 is important. 
As discussed earlier, u?o is actually a frequency resolution related parameter. Suppose 
we know the minimum distance between two frequencies is Au>. u> 0 should not be 
too much larger than Au> otherwise the resolution will not be high enough and closer 
frequencies may not be distinguished by the algorithm, u o is related to the length 
of each shot. So the above argument means that the measuring time interval should 
not be too small. Selection of w 0 should not be too small either. Too small will 
cause the algorithm to require picking a large M for covering a sufficiently large 
frequency band which will not help in increasing the frequency resolution. Based 
upon our simulation we have found that the optimal choice of aio is around the range 
of 0.4Au; ~ 2.5Ao;. 


Therefore the time interval used for each shot formulation should not be too large. 
If we get a large measuring time interval T, the best way to utilize this large data set 
is to divide it into several shots for making the length of each shot appropriate. So 
we can construct a more robust algorithm as well as obtaining the optimal resolving 
frequency u> 0 . 

For testing the performance of using different resolving frequencies, we sampled the 
signal up to T = 3 seconds. We divide this three second signal into i disjoint sections. 
i — \ means that we use whole data set as one shot: i — 2 means that we separate the 
whole data set into two shots, and etc. The resolving frequencies utilized for i disjoint 
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separations is equal to 27ri/3. We chose i = 1,2, . . . , 7 and tested the algorithm in the 
different SNR circumstances. Fig.5.4(a-c) show the simulation results for frequency 
u >2 and u> 3 . From the figures we can see that the best result is obtained by selecting u?o 
around 8.5(rad./sec.), i.e., dividing the whole data set into four disjoint shots. Also, 
the selection of ujq is not very sensitive, especially in high SNR cases. 



Fig. 5. 4(a) Variance of estimation of against resolving frequency w 0 in different 
SNR cases. 
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Fig. 5.4(b) Variance of estimation of w 2 against resolving frequency w 0 in different 
SNR cases. 




Fig. 5. 4(c) Variance of estimation of w 3 against resolving frequency w 0 in different 
SNR cases. 
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Example 2. In the second example, the received signal is composed of four 
sinusoids. The analog angular frequencies, amplitudes and phases are given by: 
u ai i = 40, u> a ,2 = 50, u> a ,3 = 60, u aA = 70; A\ = A 2 = A 3 = A 4 = 5; 4>\ = 
0, <f > 2 = 0.872, <f> 3 = 1.537, <p 4 = 1.975. The sampling frequency f 3 is 100(Hz). The 
measuring time T is 1.28 seconds. So we have total of N = 128 points. The nor- 
malized frequencies are 0.06366(Hz), 0.07958(Hz), 0.09549(Hz) , 0.11141(Hz). A is 
selected as 4. a> 0 = 2n/T = 4.91. We utilize one-shot time estimation and let m, = i 
for simplicity. There are four unknown parameters. M is chosen to be 8. Reasoning 
the same as before, we select L = 20, L\ = 25 for using the HOYW equations. It is 
very good selection for the HOYW algorithm. 

Fig.5.5(a-d) give the plots of normalized variances of estimates u>,-, i = 1,2, 3, 4 
using the modulating signal method and the HOWY equations based upon the total 
of 50 Monte-Carlo runs. Also the theoretical Cramer- Rao lower bound is presented for 
easy comparison. Same as example 1, the HOYW equations cannot identify the four 
frequencies for SNR lower than 5(dB). Fig.6 shows the power spectrum of estimated 
predictive filter for SNR=0(dB). From the graphs we can see that there are only two 
high peaks around frequencies 0.4 to 0.7. Also there are two high spurious peaks at 
other frequencies. These high spurious peaks will increase the difficulty in recognizing 
the true frequencies. 
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5.6 Concluding Remarks 


We have proposed a new algorithm for high resolution frequency estimation using the 
Fourier based sinusoidal modulating function technique. A continous autoregressive 
differential equation model was used for the analysis. The complex sinusoidal modu- 
lating signal was constructed for projecting the model into linear algebraic equations. 
Comparing with the HOYW equations, the new algorithm gave better performance, 
especially for short time observation and low SNR data. 
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Chapter 6 


Deconvolution and Parameter 
Identification for Noncausal 
Nonminimum Phase ARMA Systems 
Using Inverse Cumulants 


6.1 Introduction 

Linear time-invarant system modeling, nonminimum phase system deconvolution and 
identification play a very important role in adaptive process control, system theory, 
signal processing and data communication. Methods abound for system modeling 
and identification via the second order statistics, i.e., autocorrelation and power spec- 
trum [Sode89, Ljun87, BoxG70], and inverse autocorrelations can be used to estimate 
the coefficients of MA models driven by random Gaussian input processes [Chat79, 
Clev72]. Although the first and second order statistics which characterize a random 
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Gaussian process are quite useful tools in many pract ical applications, the problem of 
using only second order moments is that it is phase blind. Hence only the minim 
phase, power spectral equivalent modes of the system can be identified. 

Recently, some researchers have used higher order statistics ( bispectrum, trispec- 
trum, etc. ) for estimating the parameters of linear rational transfer functions [Bril67, 
Gian89, NikiST, Frie89] . It is well-known that high order statistics are phase sensitive, 
and contain information of the true phase response as well as magnitude response of 
the linear system. Also, the high order statistics are insensitive to Gaussian noise, 
which is a very important property and has many practical applications, e.g., signal 
detection in large Gaussian white noise background. Comparing to the second order 
moments, the high order statistics can be employed for identifying and reducing non- 
minimum phase systems, estimating the parameters of non-Gaussian processes, and 

modeling nonlinear systems. 

In this chapter, we focus on the issues of system deconvolution and parameter 
identification. Currently, Chiang and Nikias [Chia90, NikiST] use a noncausal AR 
process to approximate a finite FIR system and construct certain linear equations 
with respect to the parameters of an AR process based upon the third order cumu- 
lants. They invoke some adaptive gradient-type algorithms to solve these equations 
in obtaining estimates of the parameters. The problem is that their method cannot 
be applied to infinite HR systems. 

Here we are interested in the deconvolution and parameter identification of a gen- 
eral SISO noncausal nonminimum phase ARMA system driven by a random non- 
Gaussian process. Although in system theory and signal processing, causal system 
modeling is commonly used, noncausal ARMA system will also be necessary for char- 
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acterizing many processes in practice. Except for the very special case of a symmetric 
noncausal system, the high order statistics will play a crucial rule in identifying the 
system ( phase response as well as magnitude response ). 

We introduce the inverse polyspectra and inverse cumulants which are reciprocals 
of the polyspectra and cumulants of an original process. It is a natural generaliza- 
tion of the definition of inverse power spectrum and inverse autocorrelation [Chat 79, 
CIe\72]. It will be demonstrated that there is duality between cumulants and inverse 
cumulants which corresponds to the duality between the AR modes and MA modes 
of the system. For a general noncausal ARMA system, the inverse filter has been 
constructed by utilizing a noncausal AR model to approximate the original process. 
Then a relationship can be established between the inverse cumulants and the param- 
eters of the inverse filter which is unique to within a scale factor. For achieving good 
numerical performance, we use a gradient type nonlinear optimization algorithm to 

minimize an error function for matching the inverse cumulants and the parameters of 
the AR model. 

The algorithms are proposed to estimate the inverse cumulants using both fre- 
quency domain and time domain formulations. The advantage of estimating the 
inverse cumulants in the frequency domain is that the numerically efficient FFT tech- 
nique is available and no equations need to be solved. But the price we pay is that 
the algorithm is less robust in the low signal-to-noise ratios(SNR) situation. For es- 
timating the inverse cumulants in the time domain, some linear equations in terms 
of inverse cumulants have been derived for using the least squares technique. We use 
the forward regression orthogonal algorithm suggested in [Bill89] to solve the least 
squares problem because only part of the unknowns need to be estimated. The algo- 
rithm can be used to estimate a subset of the parameters, and computation time can 
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be saved. With adoption of the steepest gradient type optimization algorithm, the 
parameters of the AR filter can be identified adaptively. 

This chapter is arranged as follows. In section 6.2, we review the background 
about cumulants and polyspectra and give the problem statements. In section 6.3, 
the definitions are given of the inverse polyspectra and inverse cumulants. Section 6.4 
provides some algorithms for estimating these quantities, deconvolving and identifying 
the nonminimum phase system based on the inverse cumulants estimated both in the 
frequency and time domain. In section 6.5, some computer simulation results are 
presented. Section 6.6 is concluding remarks. 


6.2 Background and Problem Statements 


In this section, we give the problem statements and present some preliminary defini- 
tion and results related to the cumulants and polyspectra before giving the definition 
of the inverse cumulants. For a rigorous definition of the cumulants, we recommend 
[Rose85] to the reader. 

The scalar output signal {y(<)} is assumed to be a zero mean stationary , discrete 

time process, described by the following noncausal linear ARMA model: 

M+ N+ 

Y a k y{t-k)= Y hu{t-k) (6.2.1) 

k=-M~ k——N~ 

where u(t) satisfies the condition (6.2.4) specified below. 

Equation (6.2.1) can also be expressed equivalently by the impulse response ver- 
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sion: 


k=+oc 

y(t) = ^2 h(k)u(t — k) (6.2.2) 

k— — oo 

The third and fourth order cumulants contain sufficient information to characterize 
the magnitude and phase of the system. For a zero mean stationary process the 
kth order cumulants of y(t ), k = 1, 2, 3, 4, can be computed via the following equations: 

c\ = E{y(i)} = 0 
c 2 (m) = E{y(i)y(i + m)} 
c|(mi,m 2 ) = E{y(i)y(i + mi)y{i + m 2 )} 

(6.2.3) 

c^(m 1 ,m 2 ,m 3 ) = E{y{i)y(i + mi)y(i + m 2 )y{i + m 3 )} 

— c 2 (mi )c 2 (m 2 - m 3 ) - c 2 (m 2 )c 2 (m 3 - mj) 

-c 2 (m 3 )c2(mi - m 2 ) 

Suppose the system (6.2.1) satisfies the following modeling assumptions. The input 
driving signal u(t) is stationary, non-Gaussian, zero mean i.i.d. , Eu(t)u(t It) = 
cr 2 S(r), and with the kth order cumulants 7^, 0 < |7^| < 00, k > 3: 

c u k {m u m 2 , .... m k - 1) = 'ltb{rn l )b{m 2 ) . . . 6(m*_ 1) (6.2.4) 

where 6(mi) is the Kronecker delta function. u(t) has a kth order flat polyspectra, 
and is not accessible in the identification scheme. 

Also, we assume that the ARMA model (6.2.1) wfith parameters 6 = { a _ jvf — , . . . , «m+, 
6 _jv-, • • • , } is generally nonminimum phase , i.e. the system has some zeros inside 

and outside the unit circle, free of pole zero cancellations. The system may have some 
poles outside the unit circle. These poles correspond to the noncausal modes. 
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(6.2.5) 


The transfer function of the system (6.2.1) is 


A: = + oc 

H(z) = £ h(k)z- k 

k —~00 


Eg-*- b k z- k 

Eg-M- ^Z~ k 


Biz} 

A{z) 


The region of convergence is a ring including the unit circle po < \z\ < P\- po is the 
maximum amplitude of the poles inside the unit circle, p\ is the minimal amplitude 
of the poles outside of the unit circle. We assume that the system does not have any 
poles and zeros on the unit circle, i.e. A(z)B{z) / 0 for \z\ — 1 . 

The problem is to identify the system transfer function //(e^) ( phase response as 
well as magnitude response ) and reconstruct the input signal u{t) using only output 
data y(t). The method we utilize here is to introduce the inverse cumulants and 
approximate the original noncausal ARMA system via the AR model. There exists 
a direct relationship between the inverse cumulants and the parameters of the AR 
model, and the inverse cumulants can be estimated easily in the frequency domain 
using an FTT algorithm and in the time domain by solving a least squares problem 
using an efficient forward regression orthogonal algorithm. By this means, the original 
system can be characterized exactly by the AR model, and the input signal can be 
estimated by deconvolving this inverse filter with output signal y{t)- In the following, 
we present some basic facts about cumulants for later use. 


Denoting the kth order cumulants of y(t) by 4(mi,m 2 ,. . . , m^-i) , and assuming 
that E mi rnfc—i |cj(mi,m 2 ,...,m*_i)| < then the kth order polyspectra exists 
and is defined as: 


S% (u;i , U>2 , • • • ? Wk-\ ) — 


771 1 , 


+ OC 

,,mk_i = -oo 


E Jc — 1 


( 6 . 2 . 6 ) 
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Because u(t ) is a kth order white, it is easy to verify from ('6.2.4) that : 


S k (u i,a> 2 , • • • , cc'jt—i ) = Jk — constant 


(6.2.7) 


The kth order cumulants c^(mi,m 2 , . . . , m*^) can be calculated in terms of the kth 
order polyspectra: 


1 

cj((m 1 ,...,m fc _i) = 7 - - ^ — J_^S y k { 


( 27 r)* 

e J Z).=i IM,m 'dw i djj2 ■ ■ .dmk-'i 


(6.2.8) 


From equ. (6.2.5), the kth order cumulants can also be represented in terms of the 
system impulse response as : 

+oo 

cJJ(m 1 ,m 2 ,...,m/t_i) =7jf ^ h(i)h(i 4- mj) . . . h(i + m fc _i) (6.2.9) 

l = — OO 


Taking the Z-transform of equ. (6.2.9), we can get the following expression : 




( 6 . 2 . 10 ) 


So, 


SU‘ I,** *»-!) = -»? 


^(^i)-B(^ 2) • • • B(Zfc_i)B((ciZ2 • • • £fc-u) 1 ) 
A{Zi)A(z2) . . . A{Zk-\)A{(z\Z2 ■ . . 2fc_i) _1 ) 


( 6 . 2 . 11 ) 


It will be seen in the next two sections that equation (6/2.11) is important in the role 
of the inverse cumulants. 


6.3 The Inverse Cumulants 

In this section, we give the basic definition of the inverse cumulants and describe 
the relations between the inverse cumulants and the original ARMA system. Its 
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application with deconvolution and parameter identification will be presented in the 
next section. 

For the system (6.2.1) or (6.2.2), with kth order cumulants (6.2.9) and kth order 
polyspectra (6.2.10) (6.2.11) , we define the inverse cumulants and inverse polyspectra 
by simply looking for the inverse of expression (6.2.10) or (6.2.11). Thus the kth order 
inverse polyspectra SI k (z\, z 2 , • • • , z k -i) is defined as . 

S y k {z u z 2 , ■ ■ . , z k -i)S I y k (z u z 2 , . ■ ■ , Zk-i) = 1 (6.3.1) 

The coefficient of z 2 m2 ■ • • of SP k (zi, : • • • , z k - 1) will be called the kth 

order inverse cumulants at lag rr»i, m 2 , . . . , m fc _i , and will be denoted by ci y k (m u m 2 , . • . 
m k -\). ci y k (m u m 2 , . . . , m k -i) and Slftzu •••, **-i) «e related by the Z-transform. 

Noting equations (6.2.8) and (6.2.11), the following equation holds. 

aJJ(m 1 ,m 2 ,...,mfc_i) = ( 2 7 r ) fc - 1 J-* " ' J- * 7 %B k (eM , e*-*, . . . , e^-> ) 

e _J 2.=i U,m, du>iduj 2 • ■ - (6.3.2) 

Because we assume that the ARMA model (6.2.1) does not have zeros and poles on 
the unit circle , the righthand sides of (6.2.8) and (6.3.2) are integrable. The inverse 
cumulants based on equ. (6.3.2) are well-defined. 

For the ARMA model (6.2.1), by interchanging MA modes with AR modes, we 
can construct another noncausal ARMA process: 

7V+ M+ 

y 6fcU>(i — k) = y on-f(^ — &) (6.3.3) 
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where the input driving signal e(t) is an i.i.d., zero mean, non-Gaussian process. It 
has the flat kth order polyspectra. 


= Ik = \ 

Ik 


(6.3.4) 


The k th order cumulants are: 


4( m i,m 2 ,. = -^6(m 1 )<5(m 2 ) . . . 6(m k -i) 

Ik 


(6.3.5) 


The output signal w(t) is a stationary, zero mean , ergodic process. It is easy to 
verify that the kth order cumulants of w(t) are the kth order inverse cumulants of 
y{t), i.e. , 

c*(mi,m 2 , = d£(m x , m 2 , . . . , m k -i) (6.3.6) 

Also, w(t) has the kth order polyspectra equal to "' e} }u k-\ y • Ak and Bk 

are the denominator and numerator of the r.h.s. of equ.(6.2.11) . 

Model (6.3.3) shows the important dual property between the cumulants and in- 
verse cumulants. It demonstrates that reversing the signal flow through the ARMA 
system will generate the kth order cumulants which are the inverse cumulants of the 
original system. 


127 


6.4 Inverse Cumulants Based Algorithms for 
System Deconvolution and Parameter 
Identification 


In the previous section we introduced the inverse cumulants. Here we use the non- 
causal autoregressive model to approximate the original system. Also, algorithms 
are presented for estimating the inverse cumulants, identifying and deconvolving the 
system both in the frequency domain and time domain. We demonstrate more details 
in the following. 

Given the output stationary signal {y(t),t — 0, 1,2, . . .} generated from the model 

(6.2.1), we can always design a noncausal AR process for characterizing the original 

model and deconvolving the received signal y(t ): 

-1 +oo 

u(t) = /3ay(t) + ^ ~ 0 ~h Piyjt ' 0 (6.4.1) 

t — ~0O 1 = 1 


Denoting: 


+oo 


«(*) = £ ft* ! 

t= — 00 


(6.4.2) 


Assuming /3 0 = 1.0, 4>(z) is the deconvolution filter. The original model (6.2.1) will 
be rewritten as a noncausal AR filter: 


1 _ B(z) 

*(*) M z ) 


(6.4.3) 


Noting the equations (6.2.11, 6.3.1, 6.4.3), we have: 

SI y k {z u z 2 ,...,z k - X ) = ^ 2 )...$((z 1 2 2 ...^_i)- 1 ) (6.4.4) 


Eqution (6.4.4) reveals the fact that within a scale ambiguity the k th order inverse 
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cumulants of y(t) are the same as the k th order cumulants of w(t) generated by the 
deconvolution filter <$(z). 

*(z) has to be truncated for practical realization. So define: 

*(*) = E 

-M 

( P-m ? 0, fo=l) (6.4.5) 

<J>(r) will approximate the deconvolution filter. Generally speaking, the choice of 
M and N will depend upon the location of the zeros of the ARMA system (6.2.1). 
The closer the zeros to the unit circle, the larger the M and N will have to be chosen 
for a good approximation. 

It follows immediatly that the parameters 0i obey the following equation: 
c *fe(mi,m 2 ,. . . ftft+m, • • • A+m*_, 

Ik i=-M 

0 < mi < m 2 < . - • < rrik-i < (M + N) (6.4.6) 

Using the symmetric properties, we can refer to the cumulants at other lags. For 
the third order cumulants, the following symmetric condition holds: 


c 3 (m 1 ,m 2 ) = c 3 (m 2 ,m. 1 ) = c 3 (-m 1 ,m 2 - mi) 

= c 3 (m 2 - mi, -mi) = c 3 (mi - m 2 ,-m 2 ) = c 3 (-m 2 ,mi - m 2 ) 

—{M + N) < m i, m 2 < (A/ + N) (6.4.7) 

The nonredundant cumulant samples lie on the triangle described by 0 < m 3 < m 2 < 
M + N. 
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To minimize numerical sensitivity and achieve good performance in estimating the 
parameters ft of the inverse filter using equation (6.4.6), we construct the following 

objective function for matching the inverse cumulants: 

| (M+N) m2 

J(0) = Y ■■■ H e(mi,m 2 ,...,m fc _,) 2 (6.4.8) 

mk_i=0 rn \ — 0 


where, 


N-m*_ i 


Y j j 

;(mi, m 2 , • • • , rnfc_i ) = c?3(mi, m-2, . - • , ?T?.fc-i) — “T Y ft ift+”n • • • ft+ m *-i 

7* i=— M 


here, ft = 1, 0 = [ft-M, . . . J-u ft, • • • , ftv,7*] T - The character y has been omitted 
for easy notation. Here, we only use a center part of the inverse cumulants for our 
matching because a larger error will be involved for the estimates of the inverse 
cumulants close to the boundary due to the boundary effects. 


L e t = 0 for j < —M or j > N. For the case of using third order inverse 
cumulants, e.i., k = 3, the derivative of J with respect to ft {-M < l < N, l ± 0), 


is: 


2 (M+N) 

3 m 2 


|^ = -27 Y Y c ( m i> m 2) 

UfJl 7722=0 7711=0 

(/5/+mi A+ 7712 + ft -7711 ft — 771 1 +771 2 + ft —m2 ft+mj —m2 ) 


(6.4.9) 


where 7 = I/73. 

Taking the derivative of J with respect to 7, and setting the equation to zero, 7 
can easily be estimated in terms of other parameters: 


2 (M + N) 


7 = 


£ m ,i 0 a 3 (m 1 , m 2 )(E;=~-M ftft+m, ft+m,) 

2(M+N) ~ 


(6.4.10) 


By this means, we can reduce the total number of unknowns by one dimension. 


Also, the initial value guess is important for gradient based schemes to achieve 
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good results. Based on (6.4.6), very roughly we have the following equation: 


ci 3 (M + N,M + j) as 

73 


(6.4.11) 


So, it is easy to see that: 


* * C 'Zm"'n, + N ' M + }) 

J = -M, -M + 1 , . . . , 0 , . . . , N 


(6.4.12) 


do = 1 can be utilized to determine the unknown constant /. Thus, we can use the 
following form to guess the initial values although sometimes it is not quite satisfac- 
tory: 


mo) _ ci 3 (M + N,M + j) 
lj ~ ci 3 (M + N,M) 


(6.4.13) 


j = -M, —M + 1 , . . . , 0, . . . , A 7 


We comment that the problem of stability does not arise because the deconvolution 
AR filter is basically noncausal. We can always identify any unstable causal poles as 
the stable anticausal poles. The only assumption about the underlying system is that 
it has no poles or zeros on the unit circle. 

In the above development, it has been assumed that the inverse cumulants are 
already available. In practice, what we can get at hand is the received signal y(t) or 
sampled cumulants. We need to estimate the inverse cumulants using these quantities. 
In the following two subsections, we provide the algorithms to calculate the inverse 
cumulants both in the frequency and time domains. 
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6.4.1 Inverse Cumulants Estimation in the Frequency 
Domain 

Suppose we are given the measurements of the signal y(t) for 0 < t < L — 1, where L 
is the sample length. Practically, we use the sample cumulants m2, . . . , ) 

being the estimates of theoretical cumulants of the random process. The estimates 
up to the fourth order are given by the following equations: 

I L—m — 1 

C2 (m) = - ^2 y(i)y(i + m), 0 < m < L - 1 

L ;=o 

^ L — 771 2 — 1 

c 3 (mi,m 2 ) = y ^2 y(*)j/(* + m ! )y(* + m 2) 0 < mi < m 2 < I - 1 (6.4.14) 

" «=o 

c 4 (m 1 ,m 2 , m 3 ) = JeSTo "*" 1 y(*)y(* + + m 2 )y(i + m 3 ) 

-c 2 (m 1 }c 2 (m2 - m 3 ) - c 2 (m 2 )c 2 (m 3 - mi) - C2(m 3 )c 2 (mi - m 2 ) 

0 < m.i < m 2 < m 3 < L — 1 


Generally speaking, we can directly use the FFT algorithm to estimate the inverse 
cumulants by noting the equations (6. 2. 6), (6. 3.1) and (6.3.2). The basic steps are as 
follows. 

First we construct the k th order spectrum estimates: 

L — l . 5-— \ fc — 1 27 rn t m t 

M/1,/2, • • • ,/fc-i) = c l (mi,m 2 ,....m t _ 1 )e N s (6.4.15) 

771 ] , , . . , 772 ^ _ ] = — 1 

where Nj is a sufficient large positive integer of power 2 , Nj > 2L — 1, /, = 
73^, j = 1 , . . . , A: — 1, rij = 0,1,. . ., Nj - 1 . Let c fc (mi,m 2 ,...,m fc _i) = 0 if any 
variable m, is outside of the range ( — L + 1,1 — 1). Obviously, we can set up a k — 1 
dimensional FFT algorithm to calculate Sfc(/i,/ 2 Thus, we can get the 
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estimates of the inverse polyspectra: 


s7*(/„/2,...,/k-i) 


1 


(6.4.16) 


Because the system does not have any poles or zeros on the unit circle, equation 


(6.4.16) will not cause any existence problem. 


The k th order inverse cumulants then can be estimated by calculating: 


cifc(m 1 , m 2 , . - - , ) = 


N t - 1 


(A7) k 


- y; 4‘Sidh-h //.-O' 1 


fc — 1 27rn.m, 


E k — j 
»=1 


(6.4.17) 




where, <p is the phase shift depending on rq. Surely, equation (6.4.17) can be computed 
bv using a k — 1 dimensional If FT algorithm. 


A variational procedure can be employed to estimate the inverse cumulants, which 
uses the received signal y(t) directly. There is no necessity to estimate the sampled 
cumulants in advance. In this approach, we split the data into segments and estimate 
the polyspectra in each segment, then average the whole estimate to achieve numerical 
robustness. The details are as follows. 


Divide the data y(t), 0 < t < L - 1, into K segments, each containing L s samples, 
so, L = KL S . Each segment is formed as: 

y (,) (0 = y( l + (* _ 

0 < t < L, — 1, l<i<K (6.4.18) 


Suppose u; (,) (t) is a selected window function, L s is an integer of power 2. The 
Fourier transform of y (, *(t) is : 

Y (i) (fi) = £ y (i) {t)w [i) {t)e~ 3 1 < i < K (6.4.19) 


«=o 
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where, /, = n, = 0, 1, . . . , L s - 1. The FFT algorithm can be used to calculate 
the in each segment. 


Taking the average, we obtain the form of the bispectrum estimates as: 

s 3 (/„/ 2 ) = j- - h) 

A L s t=ll 


(6.4.20) 


Thus, the third order inverse cumulants can be calculated using the IFFT as: 


ci 3 (mi,m 2 ) = 


Ls — 1 

£ 


E 2 2tt n, 

t = l L s 


(L.f n^=0S 3 (fuf2) 


(6.4.21; 


For estimating the fourth order inverse cumulants, we calculate the periodogram 


for each segment: 

/'■'(/,) = f 1 r w (/,) | 2 

A3 


and computing: 


*(/,)= E sin(2t ‘ ~/ M 

" Sin7T// 


£ 

t=-L s + l 


- /2 - /a) 

As 

Taking the average for all segments: 


/(/.) = * 

A 1=1 

G(fuh,h) = A J^G (,) (/i,/ 2 ,/3) 

A ,=1 


(6.4.22) 

(6.4.23) 

(6.4.24) 

(6.4.25) 

(6.4.26) 


Then the trispectrum takes the form: 


5 4 (/i,/2,/3) = GifufiJs) - I{h)I{h)X{h + / 3 ) 

-mi{h)x{h + a ) - mnwvi + h) ( 6 - 4 . 27 ) 

Finally, the fourth order spectrum will be estimated and the fourth order inverse 
cumulants can be calculated via the IFFT algorithm. 
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All these FFT based algorithms are simple, efficient, and easy to implement. No 
linear equations need to be solved of any kind. 

6.4.2 Inverse Cumulants Estimation in the Time Domain 

Noting equation (6.3.1), and taking the inverse Z transform on both sides, we have: 

+ OO 

Y1 c k(m\ - li,m 2 - l 2 , . . . - /*_i)eijt(/i, / 2 , . . .,l k - 1 ) 

h /fc_i=-oo 

= 6(mi,m 2 ,...,m*_j) (6.4.28) 

In practice, we have to approximate the infinite summation by truncation. Based 
upon the third order inverse cumulants, let 

^ = {(ij) I ~{M + N) < i < M + TV, -(M + AO <j< M + N] 

S c = {0, j) | i < —{M + N) or i > (M + TV) or j < — (M + N) or j > (M + A r )} 
Define the truncation error: 

£(m, n) = J] c 3 (m.-i,n-j)a 3 (ij)-6{m,n) 

(«'d)€S 

= ~50 c 3 (m - *,n -j)d 3 (i,j) (6.4.29) 

(»'.i)€S c 

Because of the symmetrical property, the total unknowns of ci 3 is ( M + Af + 1 H M +^ r + 2 ) _ 
We choose the range of (m,n) as S = {(m,n) | ~{M + N)<m<n<(M + N)}, so 
we have a total (M 4- N + 1)(2M 4- 2N 4-1) equations. The estimation criterion now 
is defined as the sum of the least squares of the truncation error for (m,n) £ § : 

Er = 5Z £ 2 ( m >") (6.4.30) 

(m,n)€5 
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Define: 


R\(m, n, l) = c 3 {m - l,n) + c 3 (m,n - /) + c 3 (m + /,n + /) (6.4.31) 

R 2 {m,n,l) = c 3 (m-/,n-/) + c 3 (m + /,n) + c 3 (m,n + /) (6.4.32) 

i? 3 (m, n, A:, /) = c 3 (m — A:, n — /) + c 3 (m — l,n — k)+ 
c 3 (m -f- A:, n -f A: — /) + c 3 (m + A: — /, n + A:)-f 
c 3 (m + /,n — A: + /) + c 3 (m — A - + /, n + /) (6.4.33) 

By utilizing the symmetrical property of the inverse cumulants, we can write equa- 
tion (6.4.29) in the form: 


M+N 


((m,n) = c z {m,n)ci 3 ( 0,0)+ ^ R 1 (m,n,i)ci 3 {i,0)+ 


t=i 


M+N M+N i— 1 

S i? 2 (m,n,?)cz 3 (f,i) + ,i)ci 3 (*,j) - 

t'=l i=2 j=l 

(m, n) G 5 


(6.4.34) 


Divide all unknowns into two parts: 

. 2(M + JV) 2(M 4- A r ) 

Ca = [c* 3 (0,0),a 3 (l,0),ct 3 (l,l),...,ci ( , jj )’ 

ct 3 (M + iV, 0), a 3 (M + AT, 1), . . . ,cf 3 (M + A r , M + JV)] r ; 

2 ( M + N) x . ,2(A/ + TV) , ,, 

Ct = [ c *3( ^ h 1 > 0)? a 3( ^ 1- 1, 1), • • • , 


(6.4.35) 


c* 3 (M + TV - 1,0), . . . ,d 3 (M + TV-1,M + TV- 1)] ] 


(6.4.36) 


Now collecting all the equations for (to,ti) G 5 and constructing the standard 
norm regression form for a least squares solution, we obtain: 


E r = (*C - «c) T (^C - «c) 


(6.4.37) 
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where ( = [CftCj ] 7 and V and a c can be formed by using the equ.(6.4.34). Denoting 
the integer P by P = M + N + 1, we note that is P(2P - 1) x \P(P + 1), ( and 
a c are \P{P + 1). 


Actually we are only interested in ( a for the AR parameter matching and initial- 
ization. The length of ( b and ( a are roughly the same, i.e., only half of the total 
unknowns needs to be estimated. Also, only one element of the constant vector a c is 
1 , all others are zero. In order to take advantage of this part, the forward regressional 
orthogonal algorithm proposed in [19] has been utilized because of its simplicity and 
effectiveness. It can be used as a pruning least squares algorithm and it also can take 
the numerical advantage of special structure of the constant vector a c . 

The algorithm can be realized in the following steps: 


1. Step 1: Writing the as: 


/ <*1,1 
« 2.1 


Ol,2 


02,2 


Ol,P 2 \ 

a 2,P 2 


Vop,,i 0/p li2 ... Qp u p 2 / 

where C, = (M + N + 1)(2 M + 2 N + 1), P 2 = \(M + A' + 1)(M + N + 2), 

- [1.0,0, ... , o] r , [cr, cj] T c fi ] r . 


2. Step 2. Denoting and calculating: 


»7i,i = 1 


Aj.l ~ Oj,l ( j — 1, 2, . . . , Pi ) 


9\ 


^1,1 

Eft, A l 
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3. Step 3. Computing: 


T]t,k — _p, , 2 

^J = l A j,i 


k - 1 


— a j,A: — 9i,k^j,i 

i=l 

^l,fc 

9 k ,-^P, \2 

2_/j=i 

Jb = 2,3,...,Pa, * = 1,2, . . . , A: — 1, t/m = 1 


4. Step 4. Estimating: 

Cp 2 = 5p 2 
p 2 

6 =9i- £ *7i.*Cfc 

fc=t + l 

i = P 2 -l,P 2 -2,...,/ > a-Ta + l 

Where r a is the length of ( a - 


We make the following comment that if the matrix has full rank then the esti- 
mates of the inverse cumulants ci 3 (m,n) will satisfy the normal equation of (6.4.3/) 
and are the least squares solution. Also the procedure can be modified as an adaptive 
algorithm easily. 


6.5 Numerical Simulation Results 

In this section, some typical numerical results will be shown to demonstrate the per- 
formance of our approach. For simulating our algorithms, the independent identical 
one-sided exponential distributed input driving process with zero mean and random 
antipodal input sequences are generated for using the third order and fourth order 
inverse cumulants. In each example, 8300 random input samples are produced which 
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convolve with the true LT1 model to generate the output signals. The first 54 and last 
54 output samples are discarded to remove the transient effects of the response. Zero 
mean white i.i.d. Gaussian noise, or colored noise i/(t) are added to the output signal 
y(t) to generate the various signal-to-noise ratios. The signal-to-noisc ratio(SNR) 
is defined as 10Iog 10 ( ^|^|||j ). The colored noise was produced by passing the i.i.d. 
white Gaussian noise through an MA filter T(z) = 1 + 0.5 z~ x — 1.5z~ 2 . 

Example 1. (Noncausal ARMA(2,2) model) 

The model we simulated in our first example is the second order noncausal ARMA(2,2): 

0.3752(1 + 0.8z -1 )(l + O.lr) 

(1 -0.2725*“ 1 )(1 +0.63592) 

The system has one stable causal pole at 0.2725, one stable anticausal pole at —1.5725, 
one minimum phase zero at —0.8 and one nonminimum phase zero at —10. The 
true coefficient values (0^2, 0-\, 0o, 0i, 02, 03, 04) of the deconvolution AR filter are 
(-0.1595, 1.5948, 1.000, -1.5263, 1.2210, -0.9768,0.7815). One-sided exponential i.i.d. 
input sequences are generated with a 2 = 0.630, 73 = 1.000 for estimating the in- 
verse cumulants. Fig. 6. 1(a) shows the true third order inverse cumulants ci 3 (i,j) for 
— 12 < i, j < 12 of the model. 



Fig. 6. 1(a) The true third order inverse cumulants ct'affij) for -12 < i,j < 12. 
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For using the frequency domain formulation, we divide the output samples into 
16 segments, each containing 512 points. White i.i.d. Gaussian noise, or colored 
Gaussian noise are added to the output signal to produce the various signal-to-noise 
ratios for testing the performance of the algorithm. Equations (6.4.19), (6.4.20) and 
(6.4.21) are used to calculate the inverse cumulants via 512 points 1-D FFT and 
512 x 512 points 2-D FFT algorithms. Weight signals w(t) are selected as rectangular 
windows for simplicity. For the output signals contaminated by the i.i.d. white 
Gaussian noise with SNR=30dB, Fig.6.1(b) demonstrates the estimated third order 
inverse cumulants ci 3>F (i,j) for -12 < < 12 in the frequency domain (averaged 

over 20 Monte Carlo runs). Fig.6.1(c) presents the estimated error Di 3t p(i,j) = 
chiij) - ci 3 'F{iJ). 



Fig. 6. 1(b) The estimated third order in- 
verse cumulants C23,f(6 j) in the fre- 
quency domain via the FFT algorithms. 



Fig. 6. 1(c) The error for the estimated 
third order inverse cumulants using the 
frequency domain formulation. 


For estimating the inverse cumulants using the time domain least squares formu- 
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lation, we need to estimate the third order cumulants. We average the estimates over 
all 16 segments for reducing the numerical sensitivity, i.e., 

1 ^ 

c 3 (m,n) = — ^2c[ k) (m,n) 

16 fc=i 

1 16 i 511 — m 

= E y (fc) (0j/ (fc) (* + m)yW(r + n) 

16 A=1 512 t=0 


We choose hi — 5, N — 20. The i.i.d. white Gaussian noise is added to produce 
the SNR=30dB. Fig-6. 1 (d) gives the plot of estimated inverse cumulants ci 3jT (i,j) 
for -12 < ij < 12 (averaged over 20 Monte Carlo runs ). Fig. 6. 1(e) shows the 
estimation error = cis(i,j) — 



$ 

h 


Fig. 6 . 1 (d) The estimated third order in- 
verse cumulants 013,7(1, j) in the time do- 
main via the least square algorithm. 



Fig. 6 . 1 (e) The error for the estimated 
third order inverse cumulants in the time 
domain via the LS algorithm. 


Based upon the simulation, we have found that for a large signal -to- noise ratio case, 
there is not much difference from the estimates of the inverse cumulants calculated via 


the FFT algorithms using the frequency domain formulation or via the least squares 
algorithm using the time domain formulation. But for the small SNR situations, LS 
algorithm does a better job than the FFT algorithm. 



Fig. 6. 2(a) The true parameters of the de- 
convolution filter and the estimated pa- 
rameters using the inverse cumulants cal- 
culated via the FFT algorithms. 



Fig. 6. 2(b) The true parameters of the AR 
deconvolution filter and the estimated pa- 
rameters via the third order inverse cumu- 
lants calculated in the time domain. 


After obtaining the estimates of the inverse cumulants, we use the Fletcher an d 
Powell optimization algorithm to minimize the matching error function (6.4.8) and get 
the estimates of the AR parameters. In the situation of SNR=20djB for the additive 
colored noise, the mean and the standard deviation of the parameter estimates we got 
( from the total 20 Monte-Carlo runs ) via the third order inverse cumulants calculated 
in the frequency domain are 0- 2 , P-i, fa, fa, fa, fa) = (-0.2342 ± 0.0925,1.3679 ± 
0.3482, -1.7762 ± 0.4419, 1.0189 ± 0.3827, -0.6820 ± 0.2523,0.6794 ± 0.1836). The 
estimates we obtained via the third order inverse cumulants calculated in the time 
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domain are (0. 2 , 0-u 0u 02, 03, 0 a) = (-0.2178 ± 0.0784, 1.4623 ± 0.3215, -1.7298 ± 
0.3764, 1.0787 ± 0.3521, -0.7219 ± 0.3621,0.6932 ± 0.1832). Fig. 6 . 2 (a) and 2 (b) give 
the plots of the true deconvolution parameters and the estimated parameters. So 
the time domain least squares inverse cumulant estimation slightly outperforms the 
frequency domain FFT estimation. 

Example 2.(Noncausal ARMA(2,3) model) 

The second example we are going to simulate is the noncausal nonminimum phase 
ARMA(2,3) system: 

1.2013(1 + 0.55z)(l + 0. 62 -1 ) 

H{z) = (1 +0.37z)(l -O^-^l + O^z- 1 ] 

The system has two stable causal poles at (0.62, -0.35), one stable anticausal 
pole at —2.70, one minimum phase zero — 0.6 and one nonminimum phase zero 1.82. 
The true coefficient values of the approximate noncausal deconvolution filter are 

(/?_ 2 , 0_ u 0 O , 0 , , 02, 03, 04) = (0.1331, -0.2420, 1.0000, -0.7920, 0.2946, -0.1767, 0.1060) 

In this example, the input driving process u(t) is a binary random i.i.d. signal 
taking the values {-1,1} with the equal probability 0.5. 73 = 0.0 and 74 = 1.0. The 
skewness of u(t) is zero so we utilize the fourth order cumulants and fourth order 
inverse cumulants. Fig.6.3(a) and Fig.6.3(b) present the true fourth order inverse 
cumulant lags ct 4 (i,j, 0 ) and 014 ( 1 ,), 1 ) for —8 < i,j < 8 . 
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Fig. 6. 3(a) The true fourth order inverse 
cumulants j, 0) for —8 < i,j < 8. 



Fig. 6. 3(b) The true fourth order inverse 
cumulants ct^(f, j, 1) for —8 <i } j <8. 


The M and N have been selected as 12. For the output signals contaminated by 
the i.i.d. white Gaussian noise with SNR= 30dB, Fig.6.4(a,b) show the fourth order 
inverse cumulants ci 4iF (i, j, 0) and ci 4<F (i,j, 1) for -8 < i, j < 8 estimated in the 
frequency domain via the FFT algorithms (averaged over 20 Monte-Carlo runs). 

Fig. 6. 4(c) and Fig. 6. 4(d) demonstrate the estimation error Di 4tF {i, j, 0) = d 4 (i,j, 0) — 
«4, f (*’,;,0) and Di 4<F (i,j, 1 ) = ci 4 (i,j, 1) - ci 4tF (i,j, 1). 

Fig.6.5(a,b) plot the fourth order inverse cumulants ci 4 j(i,j, 0) and d 4 j(i,j, 1) 
calculated in the time domain using the least squares algorithm. Fig.6.5(c)(d) give 
the difference Di 4>T {t,j, 0) = c* 4 (i, j,0) - d 4J (i,j, 0) and Di 4tT (i,j , 1) = d 4 (i,j, 1) - 

Cl 4 j{l 5 J i 1 ) • 
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Fig. 6. 4(a) The estimated fourth order in- 
verse cumulants ci^ y p{i, jf, 0) for —8 < 
i,j < 8 in the frequency domain via the 
FFT algorithms. 



Fig. 6. 4(b) The estimated fourth order in- 
verse cumulants ci 4 t r(i, j, 1) for —8 < 
ij < 8 using the FFT algorithms. 



Fig. 6. 4(c) The error for the estimated Fig. 6. 4(d) The error for the estimated 

fourth order inverse cumulants c? 4 (i, j, 0) fourth order inverse cumulants 1) 

for -8 < i,j < 8 in the frequency domain for -8 < ij < 8 using the FFT algo- 

via the FFT algorithms. rithms. 
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Fig. 6.5(a) The estimated fourth order in- 
verse cumulants 0) for —8 < 

i, j < 8 in the time domain via the least- 
square algorithm. 



Fig. 6. 5(b) The estimated fourth order in- 
verse cumulants 1) f° r < 

i,j < 8 using the least square algorithms. 



Fig. 6. 5(c) The error for the estimated 
fourth order inverse cumulants 0) 

for —8 < i, j < 8 in the time domain via 
the LS algorithms. 


Fig. 6. 5(d) The error for the estimated 
fourth order inverse cumulants ci^i, j, 1) 
for —8 < t, j < 8 using the LS algorithms. 
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In the case of SNR=15dB for additive white noise, the mean and standard de- 
viation of parameter estimates of the AR filter obtained from a total of 20 Monte- 
Carlo runs via the fourth order inverse cumulants calculated in the frequency do- 
main are ((3-2, 0-i, 0U P 2 , 03, 0a) = (0.0547 ± 0.0253,-0.3243 ± 0.0863,-0.6250 ± 
0.2328, 0.3669 ±0.1429, -0.2450 ±0.0953, 0.049 ±0.022). The parameter estimates we 
got using the inverse cumulants computed via the LS algorithm are (0-2, 0-i,0i, 02, 03, 
0 4 ) = (0.2178±0.0965, -0.1353±0.0826, -0.9483±0.3544, 0.2101 ±0.0673, -0.0924 ± 
0.0556, 0.0628 ± 0.0338). Fig. 6. 6(a) and Fig.6.6(b) plot the actual AR parameters 
and the estimated parameters. 



Fig. 6. 6(a) The true parameters of the de- 
convolution AR filter and the estimated 
parameters using the fourth order inverse 
cumulants calculated via the FFT algo- 
rithms. 



Fig. 6. 6(b) The true parameters of the 
deconvolution AR model and the esti- 
mated parameters using the fourth order 
inverse cumulants computed via the LS 
algorithm. 


Also in this example, we deconvolve the output signals with the estimated AR 
filters to reconstruct the input driven process. For the i.i.d. random binary input 
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sequence { — 1, 1}, we choose zero as the threshold value for detecting the input signals. 
If the data coming out from the deconvolution filter is larger than zero, we classify it 
as 1, otherwise we classify it as -1. The probability error P(e) of the classification 
is defined as the total number of the errors in the detection of the input signals 
divided by the length of the driven input data. For 8000 input data and various SNR 
cases of a total of 40 Monte-Carlo runs, Fig.6.7 demonstrates the effectiveness of our 
deconvolution scheme. For small SNR situations, there are some differences between 
the frequency domain estimation and time domain estimation. The performance of 
the latter is slightly better than the former. Also, for additive i.i.d. white gaussian 
noise and colored gaussian noise the results will change slightly. For large SNR cases, 
all these differences vanish and both estimation procedures give the same performance. 



Fig.6.7 The probability error of recon- 
structing the input random binary antipo- 
dal signals for noncausal ARM A (2,3) sys- 
tem of example 2 versus SNR of white and 
colored Gaussian noise. 


The simulation results have shown that it is possible to identify the deconvolution 




model accurately in the presence of moderate additive noise based upon the method 
proposed in this paper. Similar to other nonminimum phase system identification 
schemes, a long length of sample data is required in order to get good estimates of 
the cumulants and the inverse cumulants. Short data records and low SNR will give 
a large estimate variance. 


6.6 Concluding Remarks 


In this chapter we have introduced the inverse cumulants and inverse polyspectra 
for parameter identification and system deconvolution. The inverse polyspectra and 
inverse cumulants are generalizations of the inverse power spectra and inverse auto- 
correlation proposed in [Chat79, Clev72]. It has been shown that there is duality 
between cumulants and inverse cumulants which relates to the AR modes and MA 
modes of the system. The basic scheme is that we use the noncausal AR filter to ap- 
proximate the original system and identify this deconvolution filter directly using the 
inverse cumulants. Unlike the approach specified in [Chia90, Niki87], which can only 
be used for finite MA systems, our method can apply to the more general noncausal 
nonminimum phase ARMA systems. There is no necessity to know the order of the 
system. The only assumption about the model is that it does not have any poles or 
zeros on the unit circle. 

Also, we have presented a procedure for estimating the inverse cumulants both in 
the frequency domain using an FFT algorithm and in the time domain using a least 
squares algorithm. It is hard to make the claim relating to the issue of consistency 
of the parameter estimation because of the inverse cumulants involved. The Monte- 
Carlo numerical simulation results demonstrate the performance of our approach. We 
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point out here that it is easy to modify our algorithm to an adaptive scheme so that 
we can track the magnitude and phase response of the system and reconstruct the 
input signals in time. 
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